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CHAPTER  I 


STATEMENT  OF  THE  PROBLEM  STUDIED 


Many  time  dependent  physical  phenomena  are  characterized  by 
nonlinear  parabolic  equations,  the  solution  of  which  is  character¬ 
ized  by  a  sharp  front,  sometimes  a  discontinuity,  propagating 
through  the  solution  domain.  Among  problems  of  this  type  are  non¬ 
linear  convective  diffusion  problems  with  dominant  convective  terms, 
or  Stefan  type  problems  such  as  the  flow  of  fluids  through  porous 
media  or  the  melting  and  freezing  of  ice.  Such  problems  are  dif¬ 
ficult  mathematically  and  numerically  because  of  the  poor  regularity 
of  the  solutions.  Moreover,  the  mathematical  theory  underlying  these 
problems  and  their  approximations  is  very  much  incomplete. 

Toward  resolving  some  of  these  issues,  a  three-year  project  was 
initiated  in  1976,  designed  to  study  not  only  the  qualitative 
features  of  solutions  of  nonlinear  problems  of  this  type,  but  also 
in  developing  numerical  schemes  for  solving  such  problems.  Since 
then  the  research  has  basically  taken  two  somewhat  distinct  directions, 
i'lrst,  a  study  of  the  use  of  variational  inequalities  as  a  means  of 
formulating  time-dependent  Stefan  problems  was  initiated.  Classes  of 
problems  considered  here  include  the  one-phase  and  two-phase  Stefan 
problems  encountered  in  porous  media  applications  and,  in  particular, 
problems  of  ablation  of  metals  and  freezing  and  thawing  of  soils. 

A  variety  of  finite  element  schemes  were  developed  and  studied  for 
these  problems,  some  of  which  proved  to  be  very  effective.  Using 
variational  inequalities  as  a  basis,  some  new  numerical  methods  were 


developed  for  two  dimensional,  two-phase  Stefan  problems  with  time 
dependent  boundary  conditions.  A  variety  of  example  problems  was 
solved,  and  a  method  was  produced  which  seems  to  be  very  effective. 
Some  of  the  results  of  this  portion  of  the  study  have  appeared  or 
will  appear  soon  In  the  literature.  The  analysis  of  the  one-phase 
Stefan  problem  was  first  completed,  including  not  only  the  identi¬ 
fication  of  effective  numerical  schemes  but  also  the  development 
of  a  priori  error  estimates  for  finite  element  approximations.  The 
studies  led  to  information  on  the  qualitative  and  quantitative 
behavior  of  the  solution  and  its  regularity,  the  behavior  of  the 
error,  and  criteria  for  the  selection  of  trial  functions  for  finite 
element  approximations.  In  these  studies,  Stefan  problems  were 
considered  without  convective  terms. 

At  the  end  of  the  first  year  of  the  project  it  became  clear  that 
to  model  realistically  certain  phenomena  characterized  by  parabolic 
equations  and  the  propagation  of  fronts,  it  would  be  more  appropriate 
to  include  convective  terms  in  the  formulation.  Indeed,  the 
presence  of  dominant  convective  terms  in  convection-diffusion  pro¬ 
cesses  is  known  to  lead  to  solutions  with  fronts  and  to  notorious 
numerical  difficulties.  Toward  resolving  some  of  these  issues,  a 
theoretical  analysis  was  initiated  to  study  the  behavior  of  highly 
non-linear  parabolic  equations  which  contained  convective  type  terms. 
Here  a  study  of  the  theory  of  evolution  problems  characterized  by 
pseudo-monotone  operators  was  performed.  Existence  theorems, 
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uniqueness  theorems,  regularity  theorems,  and  stability  results 
were  derived  for  operators  of  the  form  A(u)  +  c|u|^|Vu|r  ,  where 
q  and  r  take  on  values  appropriate  to  make  the  operator  pseudo¬ 
monotone,  and  A  is  a  non-linear  monotone  operator.  Existence 
theorems  using  methods  of  elliptic  regularization  were  also 
investigated.  Finally,  a  theory  of  Faedo-Galerkin  approximations 
and  semi-discrete  Galerkin  approximations  was  devised  and  applied 
to  finite  element  approximations  of  these  equations.  A  priori 
error  estimates  were  obtained  and  guidelines  for  the  development 
of  appropriate  numerical  methods  were  established. 

In  recent  months,  it  was  discovered  that  the  qualitative 
analysis  of  nonlinear  parbolic  problems  could  be  substantially 
generalized  to  include  the  effects  of  degenerate  coefficients  and 
to  model  such  complex  phenomena  as  two- temperature  heat  condition, 
with  degenerate  equations  and  nonlinear  convective  and  diffusion 
terms,  plus  the  effects  of  free  boundaries.  This  work  has  been 
completed  only  recently  and  required  considerable  effort.  Professors 
Oden,  Showalter,  and  Kikuchi  worked  on  this  phase  of  the  project, 
and  the  recent  work  of  Showalter  on  nonlinear  evolution  equations  has 
proved  to  be  invaluable.  We  feel  that  a  broad  theoretical  basis 
has  now  been  established  for  further  work  on  approximations  and  the 
numerical  analysis  of  this  class  of  problem. 


The  logical  extension  of  this  work  will  also  be  the  development 


of  numerical  algorithms  for  the  study  of  degenerate  nonlinear  con¬ 
vection  diffusion  problems  of  the  type  described  above  and  the 
numerical  study  of  representative  two-dimensional  problems.  Some 
encouraging  preliminary  results  have  already  been  obtained  in 
this  direction. 


CHAPTER  II 


SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS 

Important  results  were  obtained  in  four  areas: 

One :  Existence  theorems,  approximation  theorems,  a  priori 
error  estimates,  numerical  schemes,  and  finally  computer  codes 
were  developed  for  the  analysis  of  one-  and  two-dimensional,  one- 
and  two-phase  Stefan  problems  characterized  by  variational  inequal¬ 
ities. 

Two:  Existence  theorems,  uniqueness  theorems,  theorems  on  the 
stability  and  asymptotic  stability  of  solutions,  and  regularity  of 
solutions  were  developed  for  a  large  class  of  non-linear,  convective 
diffusion  problems  characterized  by  pseudo-nonotone  operators. 

Three:  A  priori  error  estimates  for  Galerkin  and  Faedo-Galerkin 
approximations  (defined,  in  general,  by  finite  element  methods)  were 
established  for  nonlinear  convection  diffusion  problems  involving 
general  pseudomonotone  operators. 

Four:  Existence  theorems  were  obtained  for  a  large  class  of 
nonlinear,  degenerate  evolution  equations  with  solutions  involving 
free  boundaries.  Applications  to  porous  media  and  two-phase  Stephan 
problems  were  completed. 
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7 .  TECHNICAL  DISCUSSIONS  (APPENDICES’) 

Appendix  A 

A  One-Phase  Multi-Dimcnsion.nl  Stefan  Problem  by 
The  Method  of  Variational  Inequalities 

1.  Introduction 

Stefan's  problem  has  been  considered  by  many  authors  since 
STEFAN  [1]  formulated  in  1891  his  mathemat ical  model  of  the  phenomena  of 
a  freezing  of  soils.  One-  and  two-phase  Stefan  problems  have  been  investi¬ 
gated  by  KAMENOMONSTSKAJA  [1],  OLEINIK  [11,  FRIEDMAN  [1],  and  others.  Vari¬ 
ous  numerical  procedures  have  been  developed  by  DOUGLAS  and  GALLIE  [1], 

JAMET  and  BONNEROT  [1],  and  others.  However,  while  most  existing  methods 
are  applicable  for  one-dimensional  problems,  not  all  are  extendable  to 
multi-dimensional  problems,  since  many  are  based  on  special  characteristics 
of  one-dimensional  case.  We  mention  here  some  typical  methods. 

(1)  After  discretization  with  respect  to  the  space  variable, 
the  n-th  time  increment  At  is  obtained  by  the  "Stefan"  condition. 


dt 


IT  grad  H  |  on  the  frozen  front 


so  that  the  following  nodal  point  becomes  frozen  by  the  condition 


„  _h_ 
At 


,  ,  _n-l  _ 

I  grad  0  ]j 


(2)  By  the  Stefan  condition,  the  location  L  of  the  frozen 
front  is  obtained  at  the  time  t  =  n  At  through  the  formula 


Ln  -  Ln-1 
At 


n-1 


[  grad  0  | 


and  then  the  domain  of  ice  is  discretized  hy  appropriate  finite  element  or 
difference  methods.  Here  £  is  the  latent  heat,  L  is  the  position  of 
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the  frozen  front,  0  is  the  temperature  field,  and  I  •  )  denotes 
the  difference  of  the  left  value  and  right  value  on  the  frozen  front. 

The  first  method  seems  to  be  possible  only  for  one  dimensional 
problems.  The  second  method  is  more  general,  but  it  becomes  difficult 
for  the  case  in  which  many  disjoint  freezing  parts  exist.  In  many 
problems,  several  frozen  fronts  may  occur  simultaneously  and  one  frozen 
front  may  grow  until  it  intersects  another.  Such  phenomena  are  difficult 
to  model  on  the  discretized  domain  at  each  time  step. 

The  formulation  introduced  by  DUVAULT  [1]  enables  us  to  re¬ 
solve  the  above  difficulties,  and  has  the  structure  of  a  strict  mathe¬ 
matical  analysis.  That  is,  after  a  special  transformation  from  the 
temperature  field  to  the  freezing  index,  the  problem  formulated  by  the 
variational  inequalities  is  a  problem  in  which  the  dependent  variable  is 
defined  on  the  whole  domain.  Moreover,  the  freezing  index  u(x)  is 
expected  to  be  continuously  differentiable  on  the  whole  domain;  that  is, 
there  is  no  discontinuity  of  grad  u  on  the  frozen  front,  while  grad  0 
is  discontinuous  there.  Thus,  the  problem  can  be  solved  within  a  fixed 
domain  without  iteration.  The  importance  of  this  formulation  is  that 
the  unfrozen  part  is  identified  with  the  portion  where  the  freezing  index 
remains  zero.  That  is,  if  we  can  obtain  numerical  values  of  the  freezing 
index,  frozen  and  unfrozen  parts  can  be  distinguisehed  by  the  value  of 
u  . 

The  special  transformation  and  the  Stefan  condition  restrict 
the  freezing  index  u  to  be  non-positive  on  the  whole  domain.  This  leads 
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us  to  the  inequation  formulation  instead  of  usual  weak  forms.  This  inequa¬ 
tion  can  be  solved  by  appropriate  optimization  techniques;  for  example, 
the  projectional  S.O.R.  method,  a  penalty  method,  etc.  The  formulation  by 
variational  inequalities  due  to  DUVAULT  is  not  only  powerful  for  its  compu¬ 
tational  aspects  but  also  well-posed  for  mathematical  and  numerical  analyt¬ 
ic  aspects,  as  shown  by  LIONS  [1],  JOHNSON  [1],  CKA  and  GLOWINSKI  [1],  and 
so  on.  These  numerical  analyses  are  well-established . 

In  this  paper,  we  describe  the  method  of  variational  inequali¬ 
ties  for  one-phase  Stefan  problems  and  give  a  computational  technique  to¬ 
gether  with  various  numerical  examples.  We  compare  the  results  of  numeri¬ 
cal  experiments  with  IIKHONOV's  exact  solution  of  a  one-dimensional  case 
[1],  confirm  our  error  estimates  for  finite  element  methods  by  numerical 
experiments,  and  analyze  some  nontrivial  two-dimensional  one-phase  Stefan 
problems. 


2.  Formulation  of  One  Phase  Stefan  Problem 

2.1.  Mathematical  Model.  Let  D  C  lRn  (n  =  1,2,3)  be  an  open  domain 
whose  subset  Q  defines  the  frozen  portion.  Let  r  be  the  boundary  of 
D  .  On  negative  temperature  g(t)  is  prescribed  as  a  function  of 

time  (a  Dirichlet  boundary  condition),  and  on  Tc  we  assume  that  the  temp¬ 
erature  8  maintains  a  value  of  zero.  The  flux  is  prescribed  on  I*_  (a 
Neumann  boundary  condition).  Tq  is  the  interface  of  ice  and  water  which 
moves  with  time  t  .  The  function  t  *  S(x)  defines  the  time  when  the 
water  x  €  D  changes  to  ice.  That  is,  S(x)  denotes  the  position  of  the 


4 


interface  I'  .  Its  inverse  relation,  x  =  S  (t)  is  given  by  L(t)  . 

Then  a  mathematical  model  of  the  one  phase  Stafan  problem  is  formvilated 
as  follows  (see  also  Figure  1). 

PROBLEM  I :  For  given  g(t)  ,  0  _<  t  <  T  ,  find  {S  ^(t)  ,  0(x,t)}  such  that 


—■  =  V  •  (kVP)  in  0  , 

O  t 

0  =  0  in  D  -  0  , 


0(x,t)  =  g(t)  on  r  , 


0(x,t)  =  0  on  P 


36 

°o  =  k  —  on  rF  , 


kV6  •  VS(x)  =  i  on  Tf 


9(x,0)  =  0  in  every  D 


(2.1) 

(2.2) 

(2.3) 

(2.4) 

(2.5) 

(2.6) 


Here  0(x,t)  is  the  temperature,  k(x)  the  thermal  diffusivity,  a  the 
constant  for  the  heat  radiation,  l  is  defined  as  £.  =  LpC  where  L  is 
the  latent  heat  per  unit  volume,  p  the  density  of  the  material,  and  C 
the  heat  capacity.  | 
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The  solution  of  the  initial  boundary  value  problem  (2.1)  - 
(2.6)  involves  two  major  difficulties.  One  is  that  the  domain  of  ice  part 
is  unknown.  Another  is  that  the  gradient  of  the  temperature  6  is  not 
continuous,  which  makes  it  difficult  to  represent  the  problem  variationally. 

2.2.  Duvaut's  Transformation-  Let  us  introduce  a  special  transformation 
of  0  into  the  freezing  index  u  by 


u (x, t ) 


/■t 

0(x,i)  dr  in  0 

■  S  (x) 


u(x,t)  =0  in  D  -  ft 
(2.7) 


following  DUVAUT  [1].  As  we  mentioned  earlier,  the  new  function  u(x,t) 
and  its  first  derivatives  Vu(x,t)  can  be  shown  to  be  continuous  on  the 
whole  domain  D  ,  while  V0(x,t)  is  discontinuous  on  fo  C  D  •  We  have 


Vu(x,t) 


«  VO  (x,t)  dT  in  SI 
JS(x) 


Vu  (x ,  t )  =0  in  D  —  SI 


(2.8) 


since  the  temperature  8  is  zero  on  To  ,  i.e.,  6(x,  S(x))  =  0  .  This 
shows,  in  fact,  that  Vu(x,  S(x))  *  0  .  Furthermore, 


V  •  (kVu) 


V  •  (kVO)  dr  -  VS (x)  •  [k(x)V0(x,  S(x)) ] 

JS(x) 


Referring  to  equation  (2.1)  and  (2.5),  we  have 
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V  •  (kVu) 


•t 

S(x) 


l 


3t 


0  dt 


l 


This  shows  that  the  field  equation  in  the  ice  part  ft  becomes 


Oil 

3t 


V  •  (k7u)  +  l  in  ft  . 


(2-9) 


And  in  the  water  domain  D  -  ft  we  have 

u  *  0  in  D  -  ft  .  (2.10) 

It  is  noteworthy  that  under  the  transformation  (2.7)  the  field  equation 
(2.1)  does  not  change  its  form  with  the  exception  of  the  force  term  Z  . 
Since  the  temperature  0  is  below  zero  degrees  centigrade  in  ft  ,  the 
heat  potential  u  is  also  less  than  zero: 

u(x,t)  <0  in  ft  .  (2.11) 

Combining  the  above  considerations,  we  have 

U  (ft"  “  7  '  *k7u)  “  l]  “  0  in  D  , 

u  0  in  D  , 

I7  -  7  •  (kVu)  -  Z  <  0  in  D  . 


(2.12) 
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ffmtrxv- 


Wc  thus  obtain  a  system  (2.12)  in  which  the  unknown  domain  0 
does  not  appear  explicitly. 

Thus,  the  initial  boundary  problem  I  is  transformed  in  the  fol¬ 
lowing  form: 


PROBLEM  II: 

For  given  g(t) 

,  0  <_  t  <_  T 

,  wi 

th 

E(0)  =  0 

.  g(t)  < 

0 

for  t  €  (0,°°)  , 

(2.13) 

find  u(x,t) 

such  that 

f  3u 

u  (7t  "  v 

' 

(kVu)  -  l 

= 

0  in  D  , 

(2.14) 

u 

< 

0  in  D  , 

(2.15) 

^  -  v 
at 

•  (kVu)  -  i 

< 

0  in  D  , 

(2.16) 

u(x,t) 

g  ( t )  on  rG 

(2.17) 

u(x,t) 

a 

0  on  rc 

(2.18) 

u(x, t) 

- 

k  •—  (x,t)  on  rp 

(2.19) 

u  (x ,  0) 

- 

0  in  D 

(2.20) 

REMARK  2.1.  Tn  the  above  formulation,  we  have  assumed  that  the  initial 
state  is  saturated  by  zero  degree  water,  i.e.,  9(x,0)  *  0  in  D  .  If  the 

initial  state  is  partially  ice  and  partially  0°C  water,  then  it  can  be  mod¬ 
ified  as  follows.  Let  0  be  the  ice  domain  in  D  at  the  initial  stage. 

o 

The  heat  potential  u  is  now  defined  by 


**»+>*** 


8 


u(x, t) 


u(x, t) 


0(x,t)  dr  in  fl  -  fi 

S(x) 


5(xtr)  dr  in 


(2.7)' 


u(x,t)  -  0  in  D  -  fi 


This  transformation,  suggested  by  FRIEDMAN  and  KINDERLEHER  [1],  reduces  the 
field  equation  to 


jju 

7u 


V  •  (kVu)  +  h 


(2.9)’ 


where  h  =  ».  in  fi  -  fi  and  'n  =  6  (x)  in  fi  ,  9  (x)  is  the  initial 

o  o  o  o 

temperature  in  fiQ  .  B 

REMARK  2.2.  In  the  above  formulation,  we  have  only  considered  the  boundary 
condition 

k  rr~  -  <t  0  on  Tr 
Bn  F 

If  r_  is  located  in  the  boundary  of  the  domain  V!  ,  defined  in  REMARK 

r  '  O 

2.1,  we  may  consider  the  boundary  condition 

k  “  -  Cl  B  +  8  on  rp  (2.4)’ 


That  is,  by  integrating  from  0  into  t  in  time,  the  condition 


9 


k 


3u 

i'n 


a  11  +  B  t  on  fj. 


(2.19)’ 


is  obtained.  The  boundary  condition  (2.19)  thus  reduced  to  the  condition 
(2.19)'.  We  note  that  if  i'p  is  located  in  the  boundary  and  fi  , 

then  (2.19)'  cannot  be  used.  8S 


3 .  Variational  Formulation 


The  corresponding  variational  formulation  for  the  problem  II, 
defined  in  the  previous  section,  will  be  derived  in  this  section.  Let 
(u ,v)  be  the  "inner  product"  on  the  domain  C  ,  i.e., 


(u,v)c 


uv  dx 
'  c 


and  let  w  =  3w/Dt  . 

Suppose  that  u  satisfies  (2.13)  -  (2.20).  Let  v  be  arbitrary 
function  such  that  v  <_  0  q.e.  in  D  x  [0,T],  v  =  0  on  Tc  ,  v  =  g  on 

I'  ,  where  T  is  a  positive  real  number  which  indicates  the  time  interval 

G 

of  the  problem.  Then,  by  integration  by  parts. 


(u,  v-u)D  +  (kVu ,  V  (v  -  u) }  D  -  (£,  v-u)D 

■  (u  -  V  •  (kVu)  -  i,  v  -  u)D  +  (kVu  •  n,  v  -  u)3J) 

where  n  «  (n^,  n,)  is  the  outward  normal  unit  vector  on  the  boundary  3D 


1 
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of  the  domain  D  .  From  the  boundary  conditions, 

(u,  v-u)D  +  (kVu,  7(v-u))d  -  (£,  v-u)Q 
=  (u  -  V  •  (kVu)  -  £,  v)D  -  (ou,  v  -  u)p 

F 

>_  -  (au  ,  v  -  u)  „ 

'  F 

Here  we  have  used  the  fact  that 

(u  -  V  •  (kVu)  -  l,  v)D  >_  0 
by  (2.12)  and  v  _<  0  a.e.  in  D  "■  (0,TJ  •  Thus, 

(u,  v-u)D+  (kVu,  7(v-u))d+  (au ,  v-u>r  >  (£,  v-u)p 

F 

a.e.  in  [0,T]  (3.1) 

is  obtained.  By  integration  of  (3.1)  in  time  [0,T]  ,  we  have  the  varia¬ 
tional  problem: 

PROBLEM  III:  Find  u  €  K  si.ch  that 

fT 

{ (u,  v  -  u)  +  (kVu ,  V(v  -  u))n  +  (nu,  v  -  u)  }  dt 
J0  D  V  *F 

>  f  (£,  v-u)  dt  (3.2) 

Jo 

for  every  v  £  K  ,  with  the  initial  condition  u(x,o)  ■  0  a.e.  in  D  , 


where 


11 


K 


{v 

v 


i 

c 


Here  the  space 


£  L2(0,T;  H1 (D)) :  v  €  L2(0,T;  L2(D>)  , 

‘  g  a.e.  on  rc  *  [0,T)  ,  v  =  0  a.e.  on 
x  [0,T]  ,  and  v  <_  0  a.e.  In  Dx  [0,T]) 

-i  2 

L“(0,T;  V)  means  that  for  every  v  €  L  (0,T;  V) 

fT  “» 

1 1  v  ||  dt  <  +  “ 

•  0  V 


(3.3) 


where  ||  •  ||  is  the  norm  of  the  space  V  .  E 

THEOREM  1.  (LIONS  [1]),  also  ICHIKAWA  [1]).  Suppose  that  mes  (V  Q)  >  0  . 
Then  there  exists  a  unique  solution  u  £  K  of  the  variational  problem 
,(3.2)  in  the  set  (3.3)  such  that 


u  £  L2(0,T;  H2(D))  A  L°°(o,T;  L2(D)) 

.  2 
u  -  V  •  (kVu)  t  L  (D  x  [0,T] )  . 


The  regularity  for  the  solution  u 
enables  us  to  consider  the  problem: 


i  .  e .  ,  u  CL  (o,T; 


l2(d>) 


PROBLEM  IV:  Find  utK  such  that 


(u,  v-u)D  +  (kVu,  V(v  -  u))D  +  (au,  v-u)r  >  (l,  v-u)p 


a.e.  in  [ 0 , T ] 


(3.4) 


for  every  v  C  K  , 
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K  =  {v  £  H  (D)  :  v  =  g  a.e.  on  I'  ,  v  =  0  a.e.  on  T  , 

c>  c 


(3.5) 


and  v  <  0  a.e.  In  D} 


B 


4.  Approximation  of  Variational  Inequality 

Since  it  is  almost  impossible  to  obtain  analytical  solutions 
of  PROBLEM  III  or  other  equivalent  forms  except  for  certain  one-dimensional 
cases,  it  is  natural  that  we  consider  approximate  methods.  Here  we  describe 
Galerkin  approximations  in  space  and  discretization  by  finite  difference 
schemes  in  time.  We  use  a  finite  element  scheme  for  spatial  approximations 
since  our  problem  may  have  an  irregular  boundary.  It  can  be  shown  that  the 
finite  element  approximation  converges  to  a  solution  of  the  given  problem. 

Let  us  consider  the  finite  element  discretization  of  PROBLEM 

IV.  Let 

B(u,v)  =  (kVu,  Vv)D  +  (au,  v) 

F 

(u,v)  =  (U.V)D  >  and  (f,v)  =  (i,v)D 


Then  (3.4)  can  be  written  by 


3u 

3T  *  V’U 


+  B(u,  v-u)  _>  (f,  v-u)  for  v  (  K 

Let  V  be  defined  by 


(4.1) 


V  -  {v  €  HX(D) :  vjp  »  g(t),  v(r  ■  0} 


(4.2) 


Then  the  admissible  set  (3.5)  can  be  written 


K 


{v  £  V:  v  <  0  a.e.  in  D) 


(4.3) 
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Mow,  in  general ,  wo  want  to  construct  finite  element  approxi¬ 
mation*;  In  II  (D)  for  sum  m  0  .  The  following  arguments  are  due  to 

ODEN  and  K1KUCHI  (1].  Let  Pj  be  a  partition  of  D  into  E-subdomains 

F, 

(finite  elements)  { )  '  such  that 

P  Q~  L 


uE  II  , 
e=l  e 


(4.4) 


C\  =  4’  for  e  4  f 


where  0  denotes  the  closure  of  0  .  Let  h  -•  dia  (£1  )  .  and  h  * 

e  e  e  e 

max  (h  )  .  We  consider  a  finite  dimensional  subspace  of  Hm(D)  consisting 
of  polynomials  of  degree  k  ,  k  >  in  >_  0  .  Let  S  be  the  finite-dimen- 

M 

ip. 

sional  subspace  of  H  (D)  corresponding  to  each  P  determined  by  the 

h 

above  polynomial  spaces.  We  construct  so  that  it  is  an  approximation 

of  Hm(D)  .  Let  =  sup  (diameters  of  all  spheres  in  .  Suppose  we 

have  the  condition:  there  is  a  constant  C  >  0  such  that 

c 


h  / p  <  C  for  every  e  ,  1  <  e  <  E 

e  e  —  o  y  —  —  ’ 


then  we  assume  there  exists  a  constant  C  >  0  ,  which  is  independent  of 
u  and  h  ,  such  that 


llU-'nhUll  s  i  C  h°  HUII  r 

"  HS  Hr 

r  _>  0  ,  0  <_  s  <_  rain  {tn,r}  ,  (4.5) 


o  *  min  {k  +  l-  s,  r-s}  , 

for  every  u  £  Hr(D)  .  Note  that  H1  (D)  -*•  HS(D)  is  a  projection  of 

Hr(D)  onto  Sj  (D)  CZ  llm(D)  .  Then  the  family  {S^}  can  provide  a  basis  of 
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the  real  Hilbert  space  for 

IL  ,  .  S,  (n)  is  everywhere  dense  in  Hm(D)  (4.6) 

Omimi  h 
—  max 

Let  (K.  }  be  a  closed  convex  subset  of  (Kj  Cl  )  which  has  the  fol¬ 

lowing  properties:  for  all  v  £  K  C V  ,  a  sequence  {v^}  In  ^  can 
constructed  as 


Vj  -*•  v  c  K  strongly  for  h  -*■  0  , 


and  the  weak  limit  u  of  the  sequence  (u^l  in  also  belongs  to  K  . 

Note  that  in  general  1  K  . 

Now  in  our  one  phase  Stephan  problem,  m  *  1  (recall  (4.2)), 
and  k  =  1  is  taken.  Thus,  we  select  piecewise  linear  polynomials 
as  a  basis  of  an  N-dimensional  space.  S^D)  .  Note  that  in  (4.5), 
r  »  k + 1  *  2  .  Thus,  K,  can  be  defined  by 


N 

{y  €  Sh:  y  =  l  Vj  ^  (Vi  £  11) 

i=l 


v 


=  0 


,  v|  <  0) 
l 


(4.7) 


where  I  ,  Z  ,  Z  denote  the  sets  of  all  nodal  points  on  Tc  ,  Tc  and 
in  D  respectively.  Note  that  v |  <_  0  implies  v^  £  0  for  all  i  , 

1  ^  i  N  .  Then  CL  K  .  For  simplicity,  we  denote 


u  tl 


N 


l 

1  =  1 


u< 


(4.8) 


where  u  €  K  ,  u  £  . 


! 


"  pvpjuuwpq 
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Discretization  In  Space.  Let  u  =  ,  v  =  v  .  Substituting  these 

into  (4.1),  we  have 


Mu  tt  <vV  +  Vt'vV  •  WV  -  0 


(4.9) 


where  M. .  is  a  mass  matrix,  K. .  a  stiffness  matrix,  f  a  force  vector 

J  J  J 

given  by 


Hij  =  (*i  *  * 


Kij  ”  B(*i  •  *]>  • 


(4.10) 


fJ  =  (f  *  V  * 


We  know  that  and  are  symmetric  by  (4.10).  Then,  interchanging 

i  and  j  , 


3u . 

Mij  IT1  (vi“ui)  +  Ki3uj(vi-u1)  -  fi(vi-ui>  i  0 


(4.11) 


Discretization  in  Time.  For  the  time  direction,  we  apply  a  finite  difference 
scheme : 


n+1  _  n 

[ij  ~A~JT~1  ^i'V  +  [°  Kij  uj1+1 

+  d-6)  Iijj  “  f1<vi“u1)  1  0  , 


(4.12) 


6 


where  denotes  j-th  point  value  of  u  =  u.^^  £  at  the  time  step 

n  ,  At  the  time  difference  between  n+1  and  n  time  step,  0  ,  0  <_  0 

1  is  the  Crank-Nicoison  coefficient.  For  0=1,  (4.12)  becomes 


(9v  vh'uh> +  B<v  vi.'V  ‘  <*•  vh'V  i  0 


where  9  denotes  the  implicit  finite  difference  operator  such  that 


9u 


n+1 


n+1  n 
Uh  -  “h 
At 


(4.13) 


Inequality  (4.12)  can  be  written  explicitly  in  terms  of  u 


as  follows: 


n+1 

i 


*  n+1,  n+1.  ;  ,  n+1. 

Ki.i  Ui  (vi‘Ui  >  --  fi(vi"ui  >  * 


(A. 14) 


where 


K.j  =  Mjj/At  +  OKjj  , 


Fij  =  Mij/At  "  Cl  -  0)  , 


fi  "  fi  ‘  Fij  Uj  • 


For  the  implicit  finite  difference  scheme  ( 0—1)  in  time  and 


the  linear  finite  element  method  in  space,  the  following  error  estimate  has 
been  established  by  JOHNSON  [1],  and  also  by  ODEN  and  KIKUCH1  [1], 


»r*Ai  i 
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THEOREM  2 .  (Error  Estimate).  Let  u  be  the  solution  of  (3.4)  and  (3.5). 
Let  be  the  solution  of  the  discrete  problem  (4.12)  and  (4.7)  at  the 

n-th  step  in  time  (0  =  1.0),  and  let 

n  n 

e  =  u  -  u, 
h 

Suppose  that  the  speed  of  propagation  of  the  frozen  front  is  of  order  tQ  . 
Suppose  that  u  €  l“’(0,T:  H2(D))  ,  u  £  L2(0,T;  H1(D))  ,  and  u  -  V  •  (kVu) 
-  I  (  L°°  (0 ,  T :  l”(D))  .  Then 

Max  I !  °n  I  i  q  +  <*  ll°ni]!||  At  i  ^(h2  +  AtM) 
n 

p  =  min  (2,  3  +  2a) 

where  a  ,  c  are  constants  independent  of  u  and  un  .  8 

h 


4,  Methods  of  Optimization 


The  discrete  problem  (4.14)  on  the  closed  convex  set  (4.7), 
defined  in  the  section  3.2,  can  be  solved  by  the  projectional  pointvise 
S.O.R.  method  as  long  as  the  matrix  ,  (4.14),  is  positive  definite, 

c.f.  CEA  and  GLOWINSKI  [1]. 


Projectional  S.O.R.  Method 

(i)  Pick  up  u"+*(0)  £  ;  for  example,  set  u 

(ii)  Suppose  that  k-th  iteration  u^+*(k)  £ 


CO)  „ 
h 


0 


is  known. 


u£+J  (k  +  0.5) 

n,i 


i-1  . 


(iii)  u£+J  (k  +  1) 


(1-u)  U  (k)  +  u  -  l  K  u"  \  (k  +  1) 
h’i  jtl  ^  h>J 

N 

l  K.  .  u"+1  (k)  +  f.  /  ic.  . 

.7,,  1 3  h,j  if  (ii) 


=  Min  (O,  U[J+J  (k  +  0.5)) 


(4.15) 


The  iteration  facLor  u  is  chosen  so  that  0  <  u>  <  2  .  Its  optimal  value 
is  decided  by  numerical  experiments,  while  the  convergence  of  the  above 
algorithm  is  obtained  for  0  <  w  <  2  ,  if  K  is  positive  definite. 

The  convergence  of  the  scheme  (4.15)  is  understood  by  the 
following  criterion: 


tolerance  = 


1  (k  +  1)  -  00 

h*1  h>1 

l  \vlh\  (k+D 

i-1  1 


(4.16) 


where  is  a  positive  small  number. 


1 
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It  is  certainly  true  that  there  are  several  other  ways  to 
solve  numerically  the  problem  of  variational  inequality  (4.14).  For 
example,  the  penalty  method,  the  Lagrange  multiplier  method,  the  fixed 
point  method,  and  so  on.  Here  we  merely  mention  the  final  forms  of  the 
above  methods  which  we  employed.  Some  numerical  results  of  a  one  dimen¬ 
sional  problem  using  several  optimization  methods  are  shown  in  the  follow¬ 
ing  section  (see  Example  5.5). 


Fixed  point  method:  Instead  of  (4 . 15) ^ »  we  have 


w-» h  ^ <) <k+i> 


The  step  (iii)  in  (4.15)  is  also  applied.  Here  p,0<p<l,isa  con¬ 
traction  factor  which  strongly  depends  on  the  problem  and  the  discretiza¬ 
tion. 


4 


3.  Numerical  Examples 


Our  numerical  scheme  obtained  in  (A. 13)  is  used  for  fixed 
mesh  of  finite  elements  and  arbitrary  time  interval  (for  the  stability 
of  the  numerical  scheme,  it  is  required  to  be  small  enough  if  0  £  6  <1/2 
is  used  in  (A. 12)),  while  other  methods,  for  example,  JAMET  and  BONNEROT 
11],  DOUGLAS  and  CALL IE  [1],  employ  variable  meshes  at  each  time  step,  or 
variable  time  interval  At  for  the  fixed  mesh  to  obtain  the  position  of 
the  frozen  front  using  the  Stefan  condition.  The  latter  method  can  be 
applied  only  for  one  dimensional  problems.  Why  we  can  use  the  fixed  mesh 
is  that  the  freezing  index  u  and  its  derivative  grad  u  is  continuous 
on  the  whole  domain  including  the  region  of  ice,  water  and  the  frozen 
front,  as  discussed  in  Section  2.  Thus  we  can  construct  the  variational 
form  and  its  approximation  without  any  restriction  of  dimension  of  space. 

5.1.  One  Dimensional  Case.  We  have  to  check  the  validity  of  the  formula¬ 
tion  of  variational  inequalities  compared  with  some  exact  solutions,  since 
for  one  dimensional  problem  analytical  solutions  are  known,  see  for  example, 
TIKHONOV  [1]. 

Suppose  that  the  following  one  dimensional  problem  is  con¬ 
sidered. 


20 

5t 

2 

K  ;  2 

Ox 

» 

x  <E  [0,1] 

e(x,o) 

=  0 

V 

x  €  [0,1] 

;  initial  condition 

(5.1) 

0(0. t) 

=  -1 

V 

t  <E  1R+ 

;  boundary  condition 

k  VC  •  VS (x) 

=  l 

on 

r 

0 

;  Stefan  condition 

22 


2 

The  thermal  diffusivity  k  is  given  by  k  =  ►:  .  Then  following  TIKHONOV 

11],  the  exact  solution  is  obtained  by 

j"  -1  +  C  4  (x/2,;/T  )  ,  for  x  £  cx/t 

0(x,t)  =  |  (5.2) 

[  0  ,  for  x  >  ctn 

where  <5>(x)  is  the  error  function,  C  the  constant  given  by  C  *  *(a/2ic)”* 
a  the  constant  which  determines  the  frozen  front  by  x  *  a/t  .  Note  that 
a  is  obtained  by  solving  a  transcendental  equation. 

The  exact  position  of  the  frozen  front  is  not  obtained  exactly, 
but  is  obtained  approximately  by  the  freezing  index. 

It  is  notable  that  the  formulation  by  variational  inequalities 
of  one  phase  problem  does  not  require  the  homogeneity  of  the  material  con¬ 
stants  k  and  2  ,  while  two  phase  problem  does,  see  KIKUCHI  and  ICHIKAWA 
[1]. 

Example  5.1.  Let  us  select  k  =  1.0  and  2.  =  1.0  V  x  £  [0,1]  .  We  use 

linear  finite  elements  with  the  raesh  size  h  =  0.1  .  The  time  Interval 

At  is  0.1  uniformly.  Then  at  time  t  =  0.2  ,  the  numerical  results  are 
compared  with  the  Tikhonov's  solution  in  Figure  5.1,  which  gives  good 
agreement.  Here  the  temperature  8^  at  j-th  nodal  point  on  n-th  time  step 
is  approximated  as  e"  =  (u“?  -  u"_1)/At  .  Since  the  linear  finite  element 
is  used,  the  gradient  Vu^  in  i-th  element  is  constant  in  each  element  and 
Is  obtained  by  (u^  -  u^_^)/h  .  Numerical  values  of  the  gradient  u” 
are  corresponding  with  the  Tikhonov's  solution  at  the  center  of  each  element 
exactly.  Q 
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Example  5.2.  In  the  previous  example,  ten  finite  elements  have  been  used 
for  the  discretization,  which  is  too  coarse  in  order  to  get  the  position 
of  the  frozen  front  properly,  while  the  temperature  and  the  gradient  of 
the  freezing  index  could  be  obtained  closely  enough  to  the  exact  values. 
Thus,  we  compute  the  same  model  with  fine  mesh  (h  ~  0.01  ,  i.e.,  100 
elements),  and  obtain  more  precise  position  of  the  frozen  front.  These 
results  are  shown  in  Figure  5.2.  El 


Example  5.3.  The  case  of  non-homogeneous  domain,  i.e.,  =  1.0  for 

x  £  [0,0.3)  and  k.,  =  100  for  x  £  (0.3,1)  ,  is  considered.  Figure  5.3 
designates  the  difference  of  the  propagation  speed  of  the  frozen  front 
between  the  homogeneous  and  the  non-homogeneous  region.  Since  k^  <<  k^  , 
the  front  propagates  more  rapidly  in  (0.3,1)  for  the  non-homogeneous 
case.  Even  though  the  values  of  k^  and  k^  are  very  different,  the 
projectional  S.O.R.  method  converges  within  almost  the  same  number  of 
iterations  as  the  case  of  uniform  material  domain.  B 


5.2.  Error  Norms.  Since  the  exact  solution  is  known  for  one  dimensional 
problem,  we  can  compute  discrete  error  norms,  and  observe  the  correspondence 
of  its  theoretical  estimate  given  in  THEOREM  2. 

Example  5.4.  The  model  problem  is  the  same  as  one  given  in  Example  5.1, 
i.e.,  with  uniform  k  =  C  =  1.0  .  In  Figure  5.4,  the  results  of  computa¬ 
tions  of  error  norms  are  shown.  Results  indicate  that  the  order  of  the 
error  in  H^-norm  is  exactly  h  which  agrees  with  the  theoretical  esti¬ 
mates  and  that  the  order  of  the  error  in  L^-norm  is  ® 


rrwnyrr  ~r  an 
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5.3.  Comparison  of  Several  Optimization  Schemes.  As  mentioned  in  the 
previous  section,  there  are  several  methods  to  solve  the  optimization 
problem  (4.14).  Here  we  compare  numerical  results  of  those  methods  for  a 
one  dimensional  steady  state  problem. 

Suppose  t lie  following  one  dimensional  steady  state  problem: 


u  K:  (u',  (v-uV)  ^  (1,  v  -  u)  V  v  K  (5.1) 

where  K  =  {v  (:  H*(0,1):  v(x)  s  0  a.e.  in  [0,1]  ,  v(0)  =  0.35)  ,  and 

we  choose  t  =  1  .  Its  exact  solution  is  f’tvon  by 


u(x) 


if  0  <  x  — 
"  "  /2 


0 


if 


x 


/2 


Example  5.5.  , The  domain  [0,1]  is  divided  into  20  finite  elements.  Using 
the  numerical  schemes  obtained  in  Section  4,  the  optimal  values  of  w  ,  X  , 
p  ,  etc.  for  each  method  are  obtained  as  follows: 


Projective  S.O.R. 

U)  = 

l.h 

Lagrange  Multiplier 

x 

0 . 04 

Fixed  Point 

P  = 

0.04 

It  is  notable  that  each  optimal  values  strongly  depend  upon  the  problem 
itself.  However,  X  and  p  may  be  chosen  by  the  following  criterion: 
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\,0  =  C  •  Max(U  /fi) 

where  C  *  0.01  '  0.05  ,  and  f  the  generalized  displacement  and 
force,  respectively,  at  a  certain  point. 

For  the  penalty  method  the  parameter  r  has  to  be  chosen 
small  enough  in  order  to  get  more  accurate  results,  but  it  depends  on  its 

discrete  element  length  h  ,  too.  The  selecting  criterion  of  c  may  be 

\ 

c  =  C  •  h 

-3  -4 

where  h  is  the  ncsh  size  (for  any  space  dimension)  and  C  ■  10  -  10 

Table  5.1  exhibits  the  comparison  of  the  results  by  these 
methods  via  the  exact  solution.  We  note  that  the  projective  S.O.R.,  the 
fixed  point  and  the  penalty  methods  arc  controlling  the  generalized  dis¬ 
placement  directly,  while  the  Lagrange  multiplier  method  is  controlling 
the  generalized  force.  According  to  the  results,  only  the  Lagrange  multi¬ 
plier  method  does  not  give  the  satisfactory  result.  F.ven  though  we  iterate 
400  times,  the  tolerance  is  bigger  than  1.0F.-5  in  the  Lagrange  multiplier 
method . 

The  projective  S.O.R.  method  is  most  effective  among  them  in 

this  problem. 

2 

The  converging  rate  0(c)  of  the.  penalty  method  in  L  -norm  as 

-2 

e  -*•  0  is  almost  1.0,  as  shown  in  Figure  5.5.  However,  for  c  <  10  , 

a  small  rate  is  obtained,  which  depends  on  the  round  off  error  of  finite 
element  discretization.  The  number  of  iterations  for  convergence  of  (4.15) 
arc  almost  the  same  (about  40  times),  that  is,  it  docs  depend  on  c  ,  if 


c  <  i  . 
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5.3.  Two  Dimensional  Case.  Since  the  formulation  and  the  numerical  pro¬ 
cedure  described  earlier  are  independent  of  the  dimension  of  the  problem, 
two  dimensional  problems  can  be  solved  without  special  considerations. 

Here  a  two  dimensional  model  is  examined  for  At  ,  w  ,  and  6  in  (4.14) 
and  (4.15).  Then  the  effect  of  lamping  of  the  mass  matrix  M  is  dis¬ 
cussed.  For  the  optimization,  the  projective  S.O.R.  method  is  employed  in 
the  following  examples  since  it  is  most  efficient  as  shown  in  Example  5.5. 

Example  5.6.  The  numerical  model  is  shown  in  Figure  5.6(a).  This  model  is 

selected  because  the  two  frozen  fronts  become  coupled  after  some  finite 

time,  so  that  any  other  methods  might  have  difficulties  to  solve  this 

problem.  Suppose  that  k  =  1.0  ,  {  =  1.0  and  h  =  1.0  ,  i.e.,  the  number 

of  finite  elements  is  10  x  10  =  100  .  The  maximum  tolerance  E  given 

c 

in  (4.16)  is  1.0F.-5  .  Let  us  fix  the  above  dimensions  in  the  following 
examples . 

For 


At  =  0.5  , 

uj  =  1.0  ,  for  the  projective  S.O.R. 

0  =  1.0  for  time  discretization 

the  numerical  results  are  shown  in  Figure  5.6(b).  B 

Example  5.7.  Figure  5.7  exhibits  the  case  At  =  0.1  at  time  t  ■  5.0  , 
which  gives  us  almost  the  same  results  as  the  case  At  **  0.5  as  shown  in 
Figure  5.6(b).  Here  8  «  1.0  and  u>  ■  1.0  are  used.  According  to  the 
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results,  the  time  increment  At  m3y  not  influence  the  numerical  results  so 
much,  if  0  =  1.0  (i.e.,  implicit  scheme)  is  chosen. 

Example  5.8.  Under  the  same  conditions  as  in  Example  5.6  except  the  relax¬ 
ation  factor  (o  ,  its  affection  to  convergence  is  checked.  We  note  that  u 
should  be  selected  in  the  range  such  as  0  •-  m  •  2  .  For  w  “  1.0  ,  1.4  , 
and  1.8  ,  the  number  of  iterations  is  given  in  Table  5.2  for  each  u  .  The 
case  of  (o  *  1.0  ,  i.e.,  the  projective  "Guass-Seidel"  method,  gives  the 
fastest  convergence.  The  calculated  temperature  field  0  is  the  same  for 
any  case  of  . 

Example  5.9.  Stability  of  the  numerical  scheme  (4.14)  is  checked  here  for 
some  Crank-Nicolson'  s  0  .  Since  the  matrices  and  are  constant 

for  each  time  step,  the  results  of  the  linear  parabolic  problem  may  be  applied; 
the  characteristic  equation  of  (4.14)  is  written  as 

det  PT[(1  +  e  At  M_1K)X  +  (-I+(l-e)At  M_1k)  )  P  =  0  (5.3) 

where  M  =  M , .  ,K=K..  ,1=6..  and  F  is  the  orthogonal  transformation 
ij  J.1  ij 

-1  T  -1 

associated  with  M  K  ,  i.e.,  P  M  K  P  reduces  to  eigenvalues  m  >  > 

. . . >  .  Then  A  in  (5.3)  is  obtained  by 

A  “  1  -  m^/(l  +  0  At  m^.) 

for  each  m^  .  In  order  that  the  scheme  (4.14)  is  stable,  it  is  required 
that  Max  )  A ^  |  <_  1  .  Clearly  A  <  1  ,  so  that  we  must  have 

(1  -  20) At  2  (5.4) 

Thus  if  0  is  selected  between  1/2  and  1,  the  scheme  is  unconditionally 
stable.  However,  for  0  <  1/2  the  scheme  may  become  unstable.  In  fact. 


0  =  0.25  and  At  =  0.5  give  an  unstable  example  as  shown  in  Figure  5.8. 


Example  5.10.  The  effect  of  lamping  of  the  mass  matrix  is  discussed. 

Let  M. .  be  a  lamped  mass  matrix  of  the  consistent  mass  matrix  M,.  defined 
ij  lj 


M.  . 
ij 


if  i  =  j 


if  i  /  j 


(5.5) 


That  is,  the  non-diagonal  terms  are  added  up  to  its  diagonal.  It  is  notable 

that  some  singular  frozen  fronts  at  the  first  time  step  t  =  0.5  are 

observed  as  shown  in  Figure  5.9(a),  which  is  considered  to  be  a  discretisa¬ 
tion  error.  However,  such  a  singular  behavior  on  the  frozen  front  is  not 
observed  if  the  lamped  mass  scheme  is  applied  as  shown  in  Figure  5.9(b). 
Furthermore,  after  several  time  steps  are  passed,  the  temperature  field 
becomes  entirely  the  same  in  both  cases  as  shown  in  Figure  5.10.  The  lamped 
mass  scheme  gives  a  kind  of  "smoothing"  effect  to  the  solution. 

The  diagonal  terms  M  of  the  lamped  mass  matrix  are  always 

greater  than  the  diagonal  terms  M  of  the  consistent  mass  matrix  (where 

1  is  not  summed) ;  the  procedure  of  lamping  always  gives  more  stability 
than  consistent  scheme  for  O'-  0  <  1/2  . 


4 
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6.  Conclusion 


We  have  shown  the  thee ry  and  appl ications  of  one  phase  Stefan 
problems  using  the  freezing  index.  Following  DUVAULT  [1],  the  problem 
described  by  the  temperature  field  has  been  transformed  to  the  variational 
inequality  in  terms  of  the  freezing  index.  Applying  finite  element  meth¬ 
ods  in  space  and  finite  difference  methods  in  time,  the  variational  ine¬ 
quality  has  been  discretized  into  a  system  of  linear  inequalities,  which 
can  be  solved  by  optimization  methods.  In  this  article,  four  kinds  of 
methods:  the  projectional  S.O.R.  method,  the  projectional  fixed  point  meth¬ 
od,  the  Lagrange  multiplier  method,  and  the  penalty  method  have  been  intro¬ 
duced  and  carefully  examined  for  their  speed  of  convergence  using  a  station¬ 
ary  problem.  Along  our  numerical  experiments,  the  projectional  S.O.R. 
method  is  the  fastest  optimization  method  among  them. 

Using  the  Tikhonov's  one  dimensional  solution,  the  numerical 
results  by  the  variatonal  inequality  have  been  compared,  and  very  close 
agreement  has  been  obtained.  Moreover,  using  the  same  one  dimensional 
example,  the  convergence  of  finite  element  methods  has  been  checked  numer¬ 
ically  for  the  case  of  linear  interpolations.  Numerical  results  have  agreed 
with  the  theoretical  one,  again. 

Several  nontrivial  two  dimensional  problems  have  been  per¬ 
formed,  and  the  choice  of  At  ,  0  ,  and  w  have  been  discussed  numerically. 
Furthermore,  the  effect  of  lumping  of  the  mass  matrix  has  been  checked. 

The  Crank-Nicolson' s  0  for  time  integration  should  be  selected  between 
1/2  and  1.  According  to  numerical  experiments,  the  iteration  factor  of  the 
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projectional  S.O.R.  method  u)  -  1.0  is  recommended.  Lumping  of  the  mass 
matrix  gives  smooth  frozen  fronts  at  the  first  few  steps  of  time 
integration. 

Thus,  the  formulation  of  one  phase  Stefan  problems  by  the  freezing 
index  is  fairly  effective  for  multi-dimensional  problems  as  shown  in  our 
discussions . 
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Figure  5.1  Comparison  of  variational  inequality  and  exact  solution  (one  dimension) 


=0.2  exact  front=0.5546 
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Figure  5.3  Front  propagation  in  non-homogeneous  domain 
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Table  5.1  Comparison  of  optimization  methods 
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Figure  5.9  (b)  Comparison  of  consistent  and  lamped  mass  matrix  at  t 


APPENDIX  B 


Numerical  Methods  for  Two-Phase  Stefan 
Problems  by  Variational  Inequalities 


Appendix  B 

Numerical  Methods  for  Two-Phase  Stefan 
Problems  by  Variational  inequalities 

1.  Introduction 

The  problem  of  freezing  end  thawing  of  two-  end  three-dimensional  ice 
fields  under  time-dependent  boundary  conditions  can  be  modelled  as  a  two-phase 
Stefan  problem,  with  the  free  boundary,  representing  the  interface  of  ice 
and  water,  unknown  a  priori.  Relatively  few  effective  numerical  methods  have 
been  proposed  for  such  problems,  and  those  which  attempt  to  treat  the  free 
boundary  by  schemes  employing  a  fixed  mesh  arc  seldom  given  a  complete  mathe¬ 
matical  justification. 

The  present  scheme  is  based  on  the  iheorv  of  Stefan  problems  using  the 
freezing  index  which  is  obtained  by  the  special  transformation  of  the  tempera¬ 
ture  field.  This  theory  is  first  introduced  by  Duvaut  [2],  and  studied  by 
Fretnond  [3].  Its  mathematical  basis  has  been  studied  by  Lions  [5]  and  Aguirre- 
Puente  and  F re'mond  [I]. 

The  purpose  of  this  article  is  to  introduce  a  numerical  scheme  for  solv¬ 
ing  two-phase  Stefan  problems  using  special  freezing  index  formulation  and 
to  discuss  its  efficiency.  While  there  are  several  mathematical  results 
availabLe  on  the  freezing  index  formulation,  such  as  theorems  on  the  existence, 
uniqueness,  and  regularity  of  solutions,  the  attempts  to  solve  it  have  been 
limited  to  one-ditnensiona  1 ;  see  Aguirre-Pucnte  and  Fremont!  [1].  We  give  here 
multidimensional  results  together  with  3ome  new  numerical  schemes. 

In  the  following  section,  the  field  equations  in  terms  of  the  freezing 
index  are  derived  1 rom  the  governing  equations  of  the  temperature  using  a 
Bpecial  trans forma tiou.  Then,  a  nonlinear  nondi f ferentiable  algebraic  system 
of  equations  are  obtained  by  discretizing  the  associated  variational  form  of 
the  freezing  index,  without  specifically  defining  the  spaces  to  which  admissible 


.  -T  ;  T-i 


1 


functions  belong  and  without  discussing  properties  of  finite  dimensional  sub¬ 
spaces  used  in  approximations.  The  nonlinear  system  of  algebraic  equations 
is  solved  by  a  modified  S.O.R.  method,  since  it  is  nondif ferentiable .  If  the 
system  is  differentiable,  the  Newton-Raphson  method  or  the  incremental  method 
may  be  applicable,  but  in  our  case  the  methods  are  not  applicable.  Moreover, 
the  nature  of  the  nonlinearity  of  the  system  implies  that  some  restrictions 
on  mesh  size,  time  increments,  and  physical  constants  such  as  the  conductiv¬ 
ities  of  the  ice  and  water  are  needed,  whereas  the  S.O.R.  method  may  converge 
without  any  such  restrictions  for  linear  systems.  We  also  discuss  some  smooth¬ 
ing  techniques  to  obtain  a  smooth  interface  oi  ice  and  water  and  give  some 
numerical  examples.  Our  numerical  experiments  indicate  that  the  freezing 
index  formulation  can  lead  to  a  powerful,  efficient,  and  simple  method  for 
solving  two-phase  Stefan  problems. 

2.  Two-phase  Stefan  Problems 

2.1  A  Mathematical  Description.  We  give  a  brief  description  of  a  formulation 
of  a  class  of  two-phase  Stefan  problems.  More  detailed  discussions  about  for¬ 
mulations  of  general  two-  (or  one-)  phase  Stefan  problems  can  be  found  in  the 
monograph  written  by  Rubinstein  [6). 

The  case  in  which  only  the  solid  or  the  melted  phase  is  governed  by 
the  heat  equation  and  the  temperature  of  the  other  phase  remains  constant,  is 
called  the  one-phase  Stefan  problem.  Tin-  two-phase  Stefan  problem  is  charac¬ 
terized  by  heat  equations  in  both  phases. 

Let  D  be  an  open  connected  subsi  t  of  I<n.  n  1.2,  3,  and  let  D  be 
divided  into  two  parts:  the  solid  part  1)^  r.nd  the  melted  part  D^.  If  the 
temperature  field  of  the  domain  at  time  t  is  f-presmted  bv  9(x,t),  then 
1)^  and  arc  defined  by 


1 


3 


D^t)  =  (x  *:  D:  6  (x,  t)  <  0 


Dg(t)  =  (x  C  D:  0(x,  t)  >0} 


The  surface  (or  maybe  subregion  of  D)  Iq  defined  by 


r0(t)  =  [x  fc  D:  0(x,t)  =  0} 


is  called  the  interface  (or  frozen  front)  of  the  solid  phase  and  the  melted 
phase  D^.  Let  the  boundary  T  of  D  be  separated  into  three  parts  and 

r3.  The  temperature  field  9 (x, t)  is  prescribed  on  the  boundary  T  .  There 
is  no  heat  flux  from  the  boundary  T^.  The  heat  flux  on  the  boundary  T,  is 
assumed  to  be  proportional  to  the  temperature  on  p  .  Then,  the  problem  can 
be  represented  by 


Cj©  -  V-  (kjW)  in  D 
C,.9  =,  V- (k  V9)  in  D„ 

C.  d 


9(x,t)  =  g(x,t)  on  P. 


3^9  (x,  t)  =  0 


on  3  g  D. 


cUHx,  t)  =  -pi9(x,  t)  +  qi(x,  t)  on  r3  HD. 


9  (x,  t)  =  0  on  rQ 

-VS  +  L  =  0  on  f. 


(2.6) 


8(x,0)  =  9q(x)  on  D 


Here  and  k^,  i  =  1,2,  are  the  mass  heat  capacity  and  the  heat  conductivity 

of  i-th  phase,  respectively,  g,  p^,  and  are  given  proper  functions,  0^ 

•  v  -  ,n 

is  the  initial  temperature  of  the  domain,  9  O0/3t,  <3.6  =.  k  ,  n  (30/dx  )  , 

’  i  i  a-1  a  a  ’ 


i  =  1,2,  n  =  (n 


l>-.->nn)  Is  the  outward  normal  unit  vector  on  T,  and 
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(2.7) 


[kve  J  =  k2(ve> ''  -  k^'Te) 


where  (^) +  is  the  limit  of  V  on  P^  coming  from  D,,,  and  (\Jr)  is  the  limit  of 


f  on  Tq  coming  from  D^.  The  function  t  r  s(x)  indicates  the  position  of 


Tp,  which  is  sometimes  written  by  x  =  L(t) .  The  value  £  is  the  latent  heat 


of  the  solid  phase. 


We  remark  that  the  portion  of  r  is  unknown  a  priori,  and  that  the 


gradient  V9  of  the  temperature  field  is  discontinuous  on  1^. 


2.2  The  Freezing  Index 

Because  of  the  discontinuity  of  the  gradient  of  the  temperature  on  the  interface 


r  ,  the  problem  cannot  be  formulated  variational ly  in  the  whole  domain  D  for 


the  temperature  field  if  the  position  of  the  interface  rQ  is  unknown.  To 


avoid  this  difficulty,  Duvaut  [2]  introduces  a  special  transformation,  which 
is;  later  called  the  freezing  index  by  Frcmond  [3): 


(2.8) 

where 


u(x,  t)  _  f  k.e  (X,  T)  d T 

6  1 


(2.9) 


i  1  if  x  £  D^C'O,  i  ~  2  if  x  (  D0(t) 


Since  k^,  i  =  1,2  are  constants,  u(x, ■)  is  differentiable  if  0  (x, •)  is  con¬ 
tinuous.  More  generally,  if  0 (x, •)  is  measurable,  u(x,-)  is  differentiable 
in  generalized  sense.  Then 


(2.10) 


1 


0(x,  t)  =  T—  U(x,  t) 
Ki 


This  implies  that 
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(2.11) 


(t)  =  (x  C.  D:  u(x,t)  <  0) 
D2(t)  (x  (.  U:  u(x,  t)  >  0 ) 


since  >0.  Under  the  assumption  that  k^  and  k,5  are  constants,  (2.3) 
implies 


V-Vu  (x,  t) 


r 


-CjQCXjO) 

-•  Cj  (x,  t) 

i  i:  X  ( 

DjCO) 

n 

i'x(t) 

-C  9(x,0) 
<  1 

+  c2(x,  t)  4  ; 

if  x  C 

Dj(0) 

n 

Dg(t) 

-CgS (X,0) 

+  c2(x,  t) 

if  x  i 

D,.(0) 

cl 

n 

Dg(t) 

^-c2e(x,o) 

4  Cx(x,  t)  -  I 

if  >:  C 

D?(0) 

n 

Dx(t) 

Applying  (2.6)  and  (2.10), 

(2.12) 

C. 

u(x,  t)  -  V-Vu(x,  t) 
i 

!l 

O 

L-». 

CD 

O 

t  e.  J 
ij 

where 

(2.13)  < 

j"  1  -  1  it  x  t  D^(0, 

[  3  s  1  if  X  C  D^O), 

i  -  2 

j  ,  2 

If  x  *  D,,(t) 

if  x  f.  Dg(0) 

(2.14)  • 

C12  =  l’  c2l  1  ~l> 

C11  * 

22  ° 

From  (2.P) , 

t 

<\j(x,  t)  /  k .  VB  (x,  r>  • 

0  1 

n(xi  d i 

t 

/  c).6  (x.  r)  d~ 

6  1 

Then,  putting 

(2.15) 

t 

g(x,  t)  =  /  k  g (x,  t)  dt, 
0  1 

1  1 

if  g(x,t)  <0,  1  =  2  If  g(x,  t)  >  0 

(2.16) 

t 

t)  ~  q  ,(x.  -Odi,  1 

0  1 

~  \  if 

X  •:  r^co,  i  =,-  2  If  X  '  d2(t), 

boundary  conditions  (2.4)  can  be  transformed  to 
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(2.17) 

u(x,  t)  = 

R(x,t) 

on  r 

(2. 18) 

c)u(x,  t) 

0 

on  I’ 

C 

and 

(2.19) 

du(x,  t) 

P. 

1 

‘  '  k. 

l 

U(X,  t)  -r  i 

whe  re 

i=l  if  x  C  Dj(t),  i  -  2  if  x  £  D.^Ct) 
Here  we  have  already  counted  the  initial  condition  of  u,  i.e., 
(2.20)  u(x,0)  0  in  D 


The  interface  conditions  (2.5)  have  been  also  taken  into  account  in  the 
above  considerations.  Therefore,  the  two-phase  Stefan  problem  (2. 3) -(2. 6) 
is  transformed  to  the  field  equation  (2.12),  the  boundary  conditions  (2.17), 
(2.18),  and  (2.19),  and  the  initial  condition  (2.20)  in  terms  of  the  freezing 
index . 


3.  Discrete  Two-phase  Stefan  Problems 

We  will  now  consider  discrete  problems  associated  with  the  problem 
((2.12),  (2.17),  (2.18),  (2.19),  (2.20))  using  finite  difference  and  element 
methods.  For  details  of  mathematical  analysis  of  the  same  class  of  two-phase 
Stefan  problems,  see  Lions  (3),  and  Agui rre-Puente  and  Fremond  [1]. 


7 


8 


3.2  Finite  Difference  Methods 

2 

Suppose  that  the  domain  D  -  K  Is  a  rectangle  which  Is  covered  by 
the  uniform  net  Zp.  Let  Z  be  the  set  of  all  nodal  points  Interior  of  the 


domain 


.  Let  Z^,  I  ,  and  Z^  be  sets  of  nil  nodal  points  rx,  rg,  and 


r  ,  respectively.  For  simplicity,  Zg  and  Z^  are  assumed  to  be  null.  Let 

the  particular  nodal  point  of  the  net  be  represented  by  the  pair  (a,©,  which 
indicates  the  position  of  the  nodal  point,  i.e.,  the  coordinates  of  the 
point  are  given  by  (OAx,  PAy)  ,  where  Ax  and  Ay  are  intervals  of  nodal 
points  in  x  and  y  directions,  respectively.  In  this  article,  for  simpli¬ 
city,  Ax  =  Ay  =  A.  Then,  the  variational  problem  (3.2)  is  reduced  to  the 
nonlinear  system: 


(3.6) 


‘  *  “ml, A  '  “a, p-1 


for  (Q£,  P)  £  Z,  and 


(3.7) 


U  rr  £ 

a,  p  fea,  p 


for  (a,p)  £  Here  the  summation  convention  is  not  applied. 

ft  a  ft  g 

Nonlinearities  can  be  included  in  and  >  but  these  make  it 

necessary  to  resolve  the  form  (3.6). 


i 
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L 


where 


*Sym)  -  k  (uS,p(m-1)  -  'O 


The  nonlinear  tenn  (n.m)i’  in  (3.9)  is  varies  like  a  step  function, 

while  other  remaining  terms  in  (3.9)  are  expected  to  change  moderately. 

This  implies  that  the  relationship 


(3.11) 


-  ~ 


.n 


it 


K,tm)  ' drF(n’ffl)iM(n)  1  1  r 


must  be  satisfied  in  order  to  get  convergence  of  the  iterative  scheme  (3.9). 

EXAMPLE  1 .  Let  us  consider  the  one-dimensional  problem,  whose  domain  D  is 
the  interval  (0,1),  material  constants  and  (i.e.,  k  ,kg,C^,  and  C^) 

are  given  constants.  Let  the  initial  temperature.  0^  be  given  by 

eQ00  ---  (x  -  1/2)  if  x  <  1/2,  -  0  if  x  >  1/2 

The  boundary  conditions  g(0,t)  and  g(l,t)  are  given  by 

8(0,  t)  *  -  1/2  g(l,  t)  t 


Then,  for  t  <  1, 


8(0,0 


-kjt/2  , 


8(1,  t)  -  k0tV2 


If  n  =  1,  and  Ci\  -  0.c>,  (3.11)  becomes 


l0-5(Vl(m)  +  lla+l  -  0.25A2d"(l,m)u^(m)  J  »0.2to2£ 

Under  the  assumption  that  u*  is  almost  the  same  witli  u^  =  9„  (i.e., 

or  a  0«  ' 

the  tine  interval  At  is  taken  to  be  sufficiently  small'*,  this  becomes 


j 


the  time  interval 


11 


0.5^  •  At  »  0.25A 2f  i.e., 

(3.12)  2k1At»A£ 

That  is,  in  order  to  satisfy  the  condition  (3.11),  the  relationships  (3.12) 
has  to  be  assumed.  Under  this  condition,  the  iterative  scheme  (3.9) 
may  converge. 

For  the  case  that  =  0.5,  Cg  -  k^  -  1.0,  i  -  100,  and  A  =  0.02, 

the  convergence  for  the  various  time  increments  At  is  obtained  in  Table  1. 

In  Figure  1,  the  position  of  the  interface  is  described  for  several 
time  increments  At.  This  shows  that  the  time  increment  At  has  to  be  small 
enough  in  order  to  treat  the  position  of  the  interface.  That  is,  it  is  pre¬ 
ferable  to  use  a  At  which  is  the  almost-limit  value  for  convergence  of  the 
iterative  scheme  (3.9). 

EXAMPLE  2.  Let  us  again  consider  a  one- dimensional, two-phase  Stefan  problem 
whose  material  constants  Cj,  Cg,  k^  kg,  and  £  arc  obtained  for  a  silty  soil 
with  twenty-. percent  moisture  content,  i.e. 

k^  =  60  k  cal/m •  day .  °C  kg  -  50  kcal/m -day. °C 

CL  =  450  kcal/ra3-°C  C?  «  600  kcal/m3-  °C 

i  =  24000  kcal/m3 

Initially,  the  soil  foundation  is  unfrozen  so  that  the  initial  temperature 

2  2 
0O  -  ax  ,  a  =  g(L,0)/l, 


0q  is  given  by 
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where  x  is  the  depth  of  the  foundation,  and  g(L,0)  is  the  initial  tempera¬ 
ture  at  the  end  depth  L.  We  specify  boundary  conditions  g(0,t)  and  g(L,  t) 
described  in  Table  2.  Let  the  depth  of  foundation  be  given  by  L  =  5m.  We 
uae  a  3-point  finite  difference  scheme  for  space  discretization  with  the  mesh 
length  A.  For  the  time  increment,  At  -  10  days  is  used. 

By  arguments  similar  to  those  in  Example  1,  the  criteria  (3.11)  be¬ 
comes 

2 

(3.13)  ||k1  §~g(0,At)  |  >>j  A2C 

Since  the  frozen  front  propagates  from  the  surface  of  the  foundation  x  =  0, 

the  condition  (3.11)  has  to  be  evaluated  at  u  1  and  n  -•  1 .  Since  (9^) 

is  almost  zero  at  x  ^  QA,  the  relation  (3.13)  is  obtained  under  the  assumption 
1  2  1 

that  is  small  enough  so  that  A  d^ijj  is  negligible.  Here  g(0,t)  is 

the  boundary  temperature  given  on  the  surface  of  the  foundation  x  0, 

We  calculate  three  cases,  i.e.,  A  0.1,  A  0.2,  and  A  =  0.5.  As 
shown  in  Table  3,  convergence  of  the  scheme  (3.9)  is  not  obtained  for  the  case 
of  A  ~  0.5.  For  A  =  0.2,  the  scheme  (3.9)  is  almost  convergent.  That  is, 
around  the  interface,  values  of  the  freezing  index  vary  periodically,  and 

-3 

the  criterion  (3.9)^^  is  not  satisfied,  choosing  c  =  10  .  However,  the 

case  of  A  =  0.1  gives  nice  stable  convergence  of  (3.9)  except  for  the  first 
few  time  steps.  For  the  first  few  steps,  the  relative  tolerance  may  not 

-3 

reach  the  given  criterion  t  ^  10  ,  since  the  value  of  the  freezing  index 

is  considerably  smaller  than  the  latent  heat  £.  Recall  that  by  a  phase  change, 
the  latent  heat  £  enters  the  force  term  in  (3.9)  .  However,  after  several 
steps,  the  global  result  becomes  stable  and  the  large  relative  tolerance  for 
the  first  few  steps  does  not  affect  the  subsequent  results.  If  we  do  not 
expect  large  relative  tolerances,  they  can  be  avoided  by  shifting  the  value 


1 
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of  the  Initial  freezing  index.  That  is,  in  place  that  the  freezing  index  is 
assumed  to  be  zero  at  the  initial  stage,  we  Just  shift  u(x,0)  to  some  posi¬ 
tive  value  Uq  whose  order  of  magnitude  may  be  the  same  as  the  latent  heat  l. 

We  also  remark  that  the  condition  (3.13)  indicates  that  the  mesh  size 
has  to  be  so  small  that  the  frozen  front  can  exceed  at  least  one  mesh  if  the 
system  is  far  from  equilibrium.  © 

We  shall  discuss  some  modifications  of  the  approximation  (3.6)  in 
order  to  get  an  efficient  method  of  solving  (3.6) . 

Ct  6  • 

First,  the  term  ^(nAt)  is  approximated  by 


(3.14) 


*%*<«*>  h 


instead  of  (3.8).  Then,  the  coefficient  d can  be  fixed  at  each  time  step. 
That  is,  material  constants  at  n-th  step  art'  determined  by  values  at  (n-l)-th 
step  which  have  been  already  obtained.  This  implies  the  modification  of 
(3.9) -(ii)  : 

<“>  •  ■£,„<">  -  (I  ?•«-«*>  * «)'  (&!,„«  ♦  vly»-i> 

J  “a.ft-lW  ‘  “«,,><!<»-»  *  ^  Vo 

2 

A  .a,0.  .  n- 1  \ 

+  aE  di  ((n_1)  At>ua,fi  j 
Na,  f3(tn)  "  c“’P(n,m)T 


(3.9)  1 


u",p(m)  =  (!■«.)«“  p(.-l)  +«^#p(-)  +f<|fp<w) 


V 


EXAMPLE  3.  Here  the  problem  described  in  Example  1  is  solved  by  the  iterative 
procedure  (3.9)'  instead  of  (3.9).  Let  Cj  0.5,  C  =  f  -  1.0,  and 

let  the  boundary  conditions  be  given  by 

g(0,t)  =  (t-l)/2,  g(l,t)  t 

The  initial  temperature  is  same  as  Example  1.  Results,  shown  in 

are  obtained  by  the  mesh  size  A  -  0.02,  the  time  increment  At  =  0.05, 

-3 

the  Iteration  factor  to  =  1.0  in  (3.9),  and  the  tolerance  e  =  10 

According  to  numerical  calculations,  see  Figure  3,  (3.9)  and  the  modified 
scheme  (3.9)  '  give  almost  the  same  propagation  of  the  solid  phase.  However, 
after  reaching  the  limit  of  propagation  of  the  solid  phase,  the 
modified  scheme  (3.9)'  gives  considerably  different  results  from  (3.9).  This 
difference  comes  from  0(d^)  -  0(£),  i.e.,  the  order  of  the  latent  heat  £ 


is  almost  the  same  as  the  one 


of  d1  -  C^/k^  and  dp  C^/kg. 


If  £  is  much 


bigger  than  d^,,  then  differences  of  results  by  (3.9)  amd  (3.9)'  may 
not  be  so  large. 

We  also  solved  the  same  problem  by  the  method  given  by  Nogi  [7 J .  Details 
of  this  comparison  are  found  in  Kikuchi  [4]. 


CC  *3 

Second,  the  term  £  of  (3.6),  which  is  related  with  the  latent 


heat  of  the  solid  phase,  is  homogenized  by 


.  ,  a,p  f  a-1,0  cul,p  a,p- 1  a, p+i  .  . 

(3.15)  Hei j  *  =  {hl<€ij  +  '  eij  +  Vj  >  +h2«ij  W<*Vh2> 

for  proper  number  hj^  and  hg.  Then,  in  (3.9)-(ii),  the  term  (n,m)  £ 

is  replaced  by  £  . 

H  lj 
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EXAMPLE  4 .  A  two-dimensional  model  is  considered  ir.  this  example.  Let 
D  =  (0,0.4)  x  (0,0.4),  C  =  0.5,  0o  =  k^  =k^  l  -  1.0.  The  initial  temperature 
is  given  by  0Q  =  0  in  D.  Boundary  conditions  are  given  by 

g(0,y)  =  (0.8-  y)  (2t  -  1) 

g(x,0.4)  (0.8  -  x)  (2t  -  1) 

g(0.4,y)  =  0.5  -/y  ( 1-3 1) 
g(x,0)  -  0.5  ,/x  (l-3t) 

The  uniform  mesh  is  employed  for  the  net  of  finite  difference  with 
A  =  Ax  =  Ay  =  0.02.  The  time  increment  At  is  0.05.  These  satisfy  the 
criteria  (3.12)  and  (3.13),  i.e., 

2kjAt  =  0.1  »  AH  =0.02 

2kjAt  =  0.1  »  A2i  =  0.0004 

At  the  9-th  and  9-th  time  step,  i.e.,  at  t  =  0.4  and  t  =  0.45,  phase  transi¬ 
tion  from  the  melted  region  to  the  solid  region  becomes  very  sensitive  in 
this  example.  At  the  7-th  time  step,  the  range  ((x,y):  0  <  x  <  0.4,  0  <  y  <  0.4, 
and  y  >  x]  is  melted  with  almost  zero  degree  temperature.  Then,  during  the 
7-th  step  to  8-th  step,  the  boundaries  (x,0)  and  (0.4,y)  become  solid.  Then, 
the  solid  phase  develops  gradually  the  boundaries  and  there  remains  the 

melted  region  inside  the  model.  We  examine  how  the  isolated  melted  region 
remains  at  t  =  0.4  and  t  =  0.45  using  the  modification  (3.15)  for  various  factors  h^, 
and  hg.  Fpur  cases  are  shown  in  Figure  4- (a)  (at  t  -  0.4)  and  Figure  4-(b)  (at  t=0.45) 
Except  for  line  of  zero  temperature,  equi-contour  lines  coincide  for 
any  choice  of  h^  and  hn  at  t.  -  0.4  and  0.45.  That  is,  the  homogenization 
(3.15)  does  not  affect  the  temperature  field  except  around  the  frozen  front. 
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However,  the  position  of  the  phase  transition  depends  strongly  upon  the 
choice  of  h^  and  h^.  Numerical  results  show  that  equi-distributed  homo¬ 
genization  (h^  s  hg  *  1)  gives  a  fairly  smooth  interface  of  the  solid  and 
melted  phases. 

Third,  the  conditions  (3.10)  may  be  replaced  by 


(3.10)  ' 


fi  =  1  if  u^3(m)  <  -c1kl 

t1  =  2  i£  >  c2k2 


cx  p 

for  the  terra  (n,m)  in  (3 . 9> .  where  e, 

numbers.  The  domain 


and  Cg  are  given  small  positive 


(3- lb)  Tr  =  (x  G  D:  -e^kj  <  u(x,  t)  <  c,,k0  } 

might  be  called  the  transient  region  of  the  solid  and  melted  phases. 


EXAMPLE  5.  The  same  example  described  in  Example  4  is  solved  by  applying 
(3.10)'  instead  of  (3.10)  in  (3,9).  Here  h^  -  0  and  h2  =  1  are  taken,  i.e.., 
no  homogenizations  are  made.  As  Example  4,  numerical  results  at  t  =.  0.4 
are  shown  in  Figure  5.  In  the  case  of  .  c,3  10  ',  shown  in  5- (c)  , 

the  whole  domain  becomes  solid.  However,  the  transient  region,  indicated 
by  small  circles,  spreads  widely.  In  this  range,  it  cannot  be  precisely  de¬ 
termined  whether  the  point  is  melted  or  frozen.  We  may  say  that  intermediate 
state  occurs  in  that  range.  The  reason  the  whole  domain  becomes  solid  is 
that  the  body  force  due  to  the  latent  heat  is  entirely  neglected  in  the  tran¬ 
sient  region.  In  the  case  of  =  c^  =  10  ^  ns  shown  in  Figure  5-(Bi,  the 

_  Q 

transient  region  certainly  becomes  narrower  than  the  case  of  =  10 
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The  temperature  fields  away  from  zero  degree  almost  coincide  in  both  cases. 
This  means  that  the  effects  of  the  modification  (3.10)*  are  limited  to  regions 
around  the  phase  change. 

3.3  Finite  Element  Methods 

For  cases  where  the  boundary  conditions  are  not  only  the  Dirichlet 
type  but  also  the  Neumann  and  the  third  types,  or  their  domains  are  con¬ 
siderably  irregular,  finite  element  discretization  is  preferable.  Let  the 
domain  D  be  triangulated.  Let  X,  be  the  set  of  all  nodal  points  in  D 

and  on  Let  iL,  and  be  sets  of  all  nodal  points  on  r  and  T  , 

2  13  13 

respectively.  Let  0^  be  the  global  interpolation  function  at  a-th  nodal 
point  which  is  constructed  by  local  interpolation  functions  (shape  functions) 
attached  to  finite  elements.  Then,  every  function  v(t)  in  H*(P)  can  be 
approximated  by 

(3.17)  v(x,  t)  =••  va(t)  $a(x) 

Here  the  summation  convention  is  applied.  In  this  section,  this  convention 
is  used  throughout--all  repeated  indices  are  summed  throughout  their  range. 

Putting  va  =  dva/dc,  (3.2)  is  discretized  by 


(3.18)  (ua)  C  Kh(t):  u'S£pvP  ,  Uas^,vP 


F^vP  i  L^vp  for  every  {vP)  t  K?(t) 

P  r3  n 


(3.19) 


Kh(t)  =  ( (va(t)  )  C  RM:  va(t)  -  |°(t) 


a  c  Ej ) 
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(3.20) 


K°(t)  {(va(t)}  C  1RN:  va(t)  0  a  i .  Z  } 


.a 


Recall  that  the  index  i  depends  upon  the  current  value  of  u  ,  and  the  index  J 


(X  q; 

depends  upon  the  initial  value  of  u  ,  i.e.,  9q.  More  precisely. 


(3.21) 


ua(t) 

<  0 

p 

l 

i  f 

0a 

eo 

<  0 

ua(t) 

>0  j 

Li- 

2 

if 

e« 

>  0 

Matrices  and  are  defined  by 


(3.22) 


>  S*  -  a,  (0.0,) 


ap  '  i'a'vC 


af'  l  cc  r 


Vectors  and  Lp  are  defined  by 


0.23)  f‘K-  <V.V' 


Here  (■,•),  a.^-,-),  and  L^(-)  have  been  defined  in  (3.1).  Then,  from  (3.18), 
it  is  necessary  to  solve  the  following  nonlinear  system 


(3.24)  . 


(ua)  €  Kh(t)  :  .  uaS^  =  +  Lj 


This  nonlinear  system  can  be  treated  by  the  iterative  algorithm  described  in 
(3.9) .  Similar  with  the  case  of  finite  difference  methods,  modifications  on 


*£p'  saf?'  and  FpJ  can  be  considered . 


First,  matrices  M  _  and  S^,  may  be  replaced  by 


(3.25) 


Mo*  (nAt)  =  (di((n-l)At)«la,0p’> 


ap 

,i 


S^OiAt)  =-  «.(C»a.0p)  ((n-l)At) 

.1  .  ..i 


Then  matrices  and  S^p  do  not  depend  upon  the  current  value  of  u  (nAt)  ,  i.e 

nonlinearity  of  matrices  are  disappeared  at  e3ch  time  step. 


19 


However,  the  homogenization  discussed  in  (3.15)  is  difficult  for 
the  finite  element  discretization.  So  we  consider  only  the  method  of 
modification  (3.10)'  instead  of  (3.10),  i.e.,  (3.21)  can  be  replaced  by 

,  i  =  1  if  ua(t)  <  -e  k 

(3.21)'  J 

i  =  2  if  ua(t)  >  e2k2 

for  the  vector  (F^). 

P 

The  numerical  scheme  which  is  employed  here  has  the  following  final 


form: 


(3.26) 


2,n 


Q!,  n 


■  - 


\  p=  1  k 


p.-ct+ 1 


Si  uP>n 
ap  k-1 


Flj)/ski 

^p  )'bca 


(a  no  summation) 


Here  u?,n  is  the  value  of  lor  the  k-th  iteration  of  the  S.O.R.  method 

k 

at  time  nAt,  u?  is  the  iteration  factor,  and 


(3.27) 


/At  +  9S 


i 

Op 


<s 


FiJ  =  FiJ 

P  P 


j  o  <  e  <  i. 


EXAMPLE  6.  Figure  6- (a)  shows  the  domain  D  and  its  boundary  conditions.  We 
employ  rectangular  linear  isoparametric  elements  with  16x  16  =  256  meshes 
whose  size  i 8  0.025,  Material  constants  are 


1.0, 

cL  =  1.0 

for 

solid  part  D^(t), 

1.0, 

Op  •-=  0.5 

for 

melted  part  D0(t) 

1.0. 
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The  initial  temperature  0^  is  given  by  L°C  everywhere.  Dirichlet  boundaries 
are  considered  at  two  points  on  the  top  surface  given  by 

f  -1.0  if  0  <  t  <0.4  and  0.5  <  t  <  0.6 
6(x,  t)  =  g(t)  -  < 

l  0.5  if  0.4'  t  <  0 . 5 

On  other  parts  of  the  boundary,  the  Neumann  condition  c6  =  0  is  assumed. 

The  time  interval  At  is  0.1.  We  use  6  =  1.0  in  (3.27),  i.e.,  the  implicit 
©-scheme  of  time  discretization  and  u'  =  1.4  as  the  overrelaxation  factor  of 
S.O.R.  method.  The  judgement  of  convergence  is  done  if  the  relative  tolerance 

Zal\’n  -  \ul'  is  loss  than  10'4- 

In  Figure  6(b)  we  show  the  case  of  =-  =  0  (see  (3. 2D  ').  Until 

time  step  3,  the  frozen  front  propagates  monotonically  with  fairly  smooth 
interface,  since  the  cooling  at  the  top  surface  is  monotone.  The  step  4 
has  somehow  unstable  values  in  temperature  C,  which  has  a  fairly  irregula  r 
shape  of  the  frozen  front.  We  think  this  is  because  the  frozen  area  is  going 
to  vanish  almost  at  this  stage.  In  order  to  avoid  this  irregularity,  we  tried 
several  cases  changing  and  e^,  which  are  shown  in  Figure  6(c) .  If  we 

compare  the  figures  at  time  step  4,  we  notice  that  the  shape  of  the  remaining 
frozen  part  is  strongly  affected  by  the  values  ot  and  e^.  However, 

fortunately,  this  frozen  area  does  not  affect  appreciably  the  temperature  field 
of  the  subsequent  time  step,  since  the  values  of  the  freezing  index  u  hardly 
change  by  the  modification  (3.10) '(the  field  variable  is  u,  not  temperature  0). 
The  selection  of  and  depends  on  our  numerical  and  experimental  experience, 
but  it  seems  to  us  that  0.001  is  iairly  proper  upon  observing 


Figure  6(c) . 
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Table  3.  Convergence  Test  for  EXAMPLE  2 


k1  ^  tg(0, Z3 1)/2 

a2  a 

convergence 

0.1 

1500. 

240. 

o.k. 

0.2 

1500. 

960. 

No  Good 

0.5 

1500. 

6000. 

No 

k1  =  60. 

J?  =  24000. 
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Figure  3.  Results  for  FXAKPI.F  3. 


Kesults  for  LXAHPLE  4  at  t=0.40 
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Figure  4-(b)  Results  for  EXAMPLE  4  at  t=0.45 
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APPENDIX  C 


Qualitative  Analysis  and  Galerkin  Approximations 
of  a  Class  of  Pseudotnonotone  Diffusion  Problems 


1.1 


1.  INTRODUCTION 

In  this  paper,  we  present  a  qualitative  analysis  of  solutions  and  Galerkin 
approximations  of  a  class  of  nonlinear  diffusion  problems  characterized  by  nonlinear 
parabolic  equations  of  the  type 

ft +  A<“>  • f 

where  u  is  an  element  of  a  separable  reflexive  Banach  space  W  densely  and 
continuously  embedded  in  another  Banach  space  V  and  A  is  a  coercive  W-pseudomonotone 
operator  from  its  domain  W  in  V  onto  the  dual  space  V’  .  Since  A  is  not 
monotone,  the  analysis  of  problems  of  this  type  is  complicated  by  the  possibility 
of  non-unique  solutions,  an  absence  of  continuous  dependence  on  the  data,  and  the 
corresponding  absence  of  stability. 

Most  of  our  attention  will  be  focused  on  the  following  class  of  problems: 

Find  u  =  u(x,t)  ,  (x, t)  £  ft  *  (0,T)  ,  such  that 

* 

+  A(u)  =  f  ,  x  £  ft  ,  0  <  t  <  T 

u  =  0  ,  x  £  9ft,  0  <  t  <  T  •  (1.1) 

u(x,0)  =  uq  , 

where 

A(u)  =  A^u)  +  A2(u) 

=  -V  •a(x,Vu)  +  b(x,u,Vu)  ,  x  £  £2  (1.2) 

a(x,Vu)  =  a(x)Vu  +  k(x)|Vu|P  ^Vu 
a,k  6  L°°(ft)  ,  2  £  p  <  °° 

a(x)  _>  aQ  >  0  ,  k(x)  >  kQ  >  0  a.e.  on  ft 

and  b(£,0  =  b(x,£,0  is  a  totally  Frechet  differentiable  function  in  IR  x  Rn 
for  which  there  are  positive  reals  q  and  r  such  that 

|b(r„0  )  1  c|r,|q|c|r 

|a5b(c,p|  <  cq | c! q_1  U I r  (q  i  o)  *  (1.4) 

i^bu.pi  <  crkiqur_1 


(r  i  0) 


1.2 


and 


q  =  0  or  q  >  1  ,  r  =  0  or  r>l 


(1.5) 


l<_q+r<p-l  | 

Here  ft  is  an  open  bounded  domain  in  Rn  with  boundary  8ft  ,  0  <  T  <  00  ,  f 
is  given  data  in  ft  x  (0,T)  ,  and  uq  is  initial  data  given  on  ft  . 

The  operator  A^  is  a  generalization  of  a  part  of  a  nonlinear  heat  diffusion 
operator  proposed  by  Coleman  and  Mizel  [4].  With  a  =  0  ,  it  also  corresponds, 
when  spaces  of  vector  functions  with  zero  divergence  are  considered,  to  the  diffusion 
part  of  the  generalized  Navier-Stokes  equation  studied  by  Lions  in  [10]  and  [11]. 

The  general  operator  may  model  a  convective  part  of  the  process.  For  example, 

convective  terms  such  as  those  in  the  Navier-Stokes  equations,  or  in  heat  conduction  prob¬ 
lems  when  convective  terms  due  to  chemical  reactions  are  present.  This  type  of  diffusion 
equation  is  important  beyond  the  study  of  thermomechanical  phenomena;  it  also  occurs 
in  modeling  biological,  social  and  other  phenomena  (Cf.  Fitzgibbon  and  Walker  [6] 
for  a  survey  of  nonlinear  diffusion  models.)  We  should  mention  Tsutsumi's  study, 

[16],  where  a  non-monotone  parabolic  problem  in  which  an  equation  of  the  form 
n  __ 

3u/3t  -  Y,  ,  3(|3u/3x.lP  3u/3x.)/3x.  -  uq  =  0  is  analyzed,  which  is  a  special 

t,i=l  l '  li 

case  of  the  problem  considered  here. 

This  study  is  divided  into  two  principal  parts.  Part  I,  Pseudomonotone 
Parabolic  Problems,  is  devoted  to  the  study  of  the  existence  of  solutions  of  a  general 
class  of  non-monotone,  nonlinear  parabolic  equations.  After  some  preliminaries 
on  properties  of  certain  function  spaces  are  laid  down  in  the  section  following 
this  Introduction,  a  general  class  of  nonlinear  parabolic  problems  involving 
coercive  pseudomonotone  operators  is  given  in  Section  3  together  with  an  existence 
theorem  for  such  problems.  An  existence  theorem  for  problems  of  this  type  has 
been  given  by  Lions  [11],  but  to  lay  groundwork  for  our  study  of  Galerkin  approximations 
taken  up  in  Section  9  of  the  paper,  we  give  an  alternate  proof.  The  principal  tool, 
both  in  the  analysis  of  the  general  problem  and  in  the  approximation  of  the  model 
problem  (1.1)  -  (1.5),  is  the  construction  of  an  elliptic  regularization  of  the 


given  problem.  This  regularization,  o£  the  type  used  by  Lions  [10],  is  introduced 
and  analyzed  in  Section  4  and  the  proof  of  the  existence  theorem  for  the  pseudomonotone 
parabolic  problem  is  completed  in  Section  5.  In  Section  6  of  the  paper,  we  give 
a  brief  summary  of  some  results  of  Oden  [13]  which  provide  sufficient  conditions 
for  pseudomonotonicity  of  nonlinear  operators  on  reflexive  Banach  spaces. 

Part  II  of  this  study  is  concerned  with  the  specific  class  of  nonlinear  diffusion 
problems  characterized  by  (1.1)  -  (1.5).  In  Section  7,  we  show  that  the  operator  A 
defined  in  (1.2)  is  coercive  and  pseudomonotone  on  a  dense  continuously  embedded 
subspace  of  the  Banach  space  LP(0,T;  W^,P(Q))  and  that  solutions  to  (1.1)  do 
exist  in  L  (0,T;  L^(fi))H  LP(0,T;  W^,P(Q))  under  the  conditions  (1.3)  -(1.5). 

In  general,  multiple  solutions  will  exist  to  (1.1)  and  there  cannot  exist  a  continuous 
dependence  on  the  data.  However,  regularity  conditions  on  the  solutions  can  be 
given  which  will  guarantee  their  uniqueness,  and  these  are  discussed  in  Section  8. 

Sections  9  and  10  are  devoted  to  studies  of  Galerkin  approximations  of  the 

model  problem  (1.1)  -  (1.5).  In  Section  9,  we  describe  properties  of  space  -  time 

Galerkin  approximations  of  solutions  in  LP(0,T;  W^,P(fi))  and  we  give  an  approximation 

theorem  which  establishes  their  strong  convergence  (in  LP(0,T;  W^,P(ft))).  We 

o 

also  derive  error  estimates  for  such  approximations.  Finally,  in  Section  10,  we 
describe  Faedo-Galerkin  (semi-discrete)  approximations.  We  note  that,  in  general, 
this  type  of  semi-discrete  approximation  is  not  necessarily  well-defined  for 
coercive  pseudomonotone  parabolic  problems.  However,  in  the  case  of  our  model 
problem,  it  is  proved  in  Theorem  10.1  that  sufficient  conditions  are  satisfied 
which  guarantee  existence  and  also  uniqueness  of  Faedo-Galerkin  approximations 
to  problem  (1.1)  -  (1.5).  We  also  prove  sufficient  conditions  for  weak  and  strong 
convergence  of  such  approximations  and  we  establish  corresponding  approximation 


error  estimates. 


2. 


PART  I.  PSEUDOMONOTONE  PARABOLIC  PROBLEMS 
2.  Some  Preliminaries 

The  following  notations  and  conventions  will  be  in  force  throughout  this 
s  tudy : 

(V,  |j  *  ||  )  =  a  real,  separable,  reflexive  Banach  space. 

(V* ,  ||  •  Ij^)  =  the  dual  space  of  V  . 

<C*,V>  ~  duality  pairing  on  V'  x  V  ;  i.e.,  for  v' €  V’  and 

v  €  V  ,  <v',v>=  v'  (v)  . 

(H,  (• , •) , | • | )  =  a  real  Hilbert  space  identified  with  its  dual,  in  which  V 
is  densely  and  continuously  embedded:  V  G|H  =  H' .  Then  H  is  a  pivot  space  such 


that 

VQIIQV'  (2.1) 

Next,  denoting  by  t  £  [0,T]  ,  0  <  T  <  °°  ,  the  time  variable,  we  introduce 
the  space  of  vector  functions  of  time 

(l/.  Ill -III)  =  LP(0,T;V) 


v:  [0,T]  •+  V  ;  |||v| 


fT 


v(t)  '  dt 


■jl/p 


<  00  }  2  <  p  < 


(2.2) 


which  is  a  separable,  reflexive  Banach  space,  whose  dual  space  can  be  identified 
as  , 

(l/'.IIHII*)  *  LP  (0,T;V)  ,  p'  =  p/(p-l)  . 

[*,•]  =  duality  pairing  on  (/’  *  1/  ,  i.e.,  for  v'  £  I/'  and 

v  e  V  , 
fT 

[v',v]  =  <v'  (t)  ,v(t)>  dt  (2.3) 

J0 

2 

(H,  (•, *)^, | • |^)  =  L  (0,T;H)  equipped  with  the  natural  inner  product  and  norm, 
which  being  identical  with  its  dual,  is  a  pivot  Hilbert  space  such  that 

1/QtfQl/'  (2.4) 

P((0,T))  =  space  of  test  functions  defined  on  (0,T)  ;  i.e., 


•321. 


i 


2.2 


<t>  €  P((0,T)) ■<*=->■  4>  €  Co((0,T))  with  the  usual  locally  convex  linear  topological 
structure. 


D'  ((0,T);X)  =  L(P((0,T)) ,X)  =  the  space  of  distributions  on  P((0,T))  with 


values  in  some  normed  linear  space  X  . 

We  observe  that  since  l/CP'((0,T);V)  ,  every 

on  fl((0,T)),  also  denoted  by  V  ,  with  values  in  v 
fT 


v(4>)  = 


v(t)(j>(t)dt  , 


V 


4>£P((0,t)) 


J0 


v  £  V  defines  a  distribution 
,  given  by 


whose  distributional  time  derivatives,  also  belonging  to  O' ((0,T);V)  ,  are  defined 
by  i-T 


4  CM  -  C-l)m 

1 


0 


v(t)  ^-1^.1  dt  ,  V  <0  6  P((0,T)) 
. .  m 


dt 


Finally,  we  introduce  the  separable,  reflexive  Banach  spaces  (Ci,  j||  •  Illy)  a°d 


(W»  III"  lllw)  defined  by 


U  =  {v:  v  6  V  ,  v  =  3v/3t  6  H} 

IIMIIU  "  Ill-Ill  +  |v|„ 

W  =  {v:  v  6  1/  ,  v  =  3v/3t  6  V' } 

ll|v||!w  -  III  v||l  +  |||v|||* 


(2.5) 


(2.6) 


which  satisfy  the  relation  (with  dense  inclusions  and  continuous  injections) 


liQUIQI/  (2.7) 

and  whose  elements  possess  the  following  properties  (cf.  [8]  and  [10]): 

i)  W  is  continuously  embedded  in  C([0,TJ;H)  ;  i.e.,  if  v6  W  ,  then, 
after  an  eventual  modification  on  a  set  of  measure  zero  in  [0,T]  ,  v 
is  continuous  from  [0,T]  into  H  and  there  is  a  constant  K  , 
independent  of  v  ,  such  that 

sup  |v(t)|  <  k|||v||L  (2.8) 

teto.T] 


ii)  If  u,v  £  W  ,  then  u,v  satisfy  the  Green's  formula 
[u,v]  =  (u(T),v(T))  -  (u(0) ,v(0) )  -  [v,u] 


(2.9) 


2.3 


iii)  The  trace  mappings  v  v(0)  and  v  h-  v(T)  from  W  -*•  H  are  surjective 
and  their  restrictions  to  U  have  range  dense  in  H  ;  i.e., 

{v(0) :  v  6  W}  =  H  =  {v(T) :  v  €  W}  (2.10) 

{v(0):  v  6  11}  and  {v(T):  v  £  U)  are  dense  in  H  (2.11) 

We  also  remark  that  we  make  frequent  use  of  Young's  inequality  in  subsequent 
analyses:  If  x,y  £  R  ,  1  <  s  <  °°  ,  s'  =  s/(s-l)  ,  and  b  is  any  real  >  0, 

then 


„  b  i  ,s  ,  1 

xy  1  —  lxl  +  — rf 


s'b‘ 


(2.12) 


3.  Existence  Theorem 


With  the  conventions  of  the  previous  section  in  force,  we  now  consider  a 

general  class  of  non-monotone  evolution  problems  characterized  as  follows: 

Given  f  €  l/'  and  u  €  H  ,  find  u  (111  such  that 

o 


f~  +  A(u)  =  f 


u(0)  =  u 


(3.1) 


Here  A  is  an  operator  from  1/  into  (/'  ,  possibly  depending  upon  the  time 

parameter  t  6  (0,T)  ,  satisfying  the  following  conditions: 

£ 

AI.  A:  (/-*■(/*  is  W-pseudomonotone  ,  i.e., 

i)  A  is  bounded  in  the  sense  that  it  maps  bounded  sets  in  V  into 
bounded  sets  in  (/'  . 

ii)  If  (u^}  Clll  is  a  sequence  converging  weakly  to  u  6  W  and  if 

lim  sup  [A(u  )  ,  u  -u]  <  0 
n  n  — 

n-*» 

then,  V  v  €  V  , 

lim  inf  [A(u  )  ,  u  -v]  _>  [A(u)  ,  u-v] 
n-*» 


All.  A:  (/-*-(/’  is  coercive,  i.e., 

^  +  00  as  IK v | 


We  will  refer  to  A:  (/  -*■  U'  as  pseudomonotone  if  it  is  U-pseudomonotone. 


3.2 


Theorem  3.1.  Let  the  operator  A:  V  V'  satisfy  conditions  (AI)  and  (All). 
Then  there  exists  at  least  one  solution  u  6  111  to  problem  (3.1).  Q 

A  proof  of  this  theorem  was  given  by  Lions  [11,  Chap.  3]  which  makes  use  of 
the  method  of  elliptic  regularization.  Lions'  method,  which  was  introduced  in  his 
study  of  linear  parabolic  problems  [9],  effectively  involves  converting  (3.1)  into 
an  elliptic  problem,  depending  on  a  real  parameter  e  >  0  ,  and  constructing 
solutions  to  (3.1)  as  limiting  cases  when  e  ■+■  0+  .  When  A  is  a  spatial 
differential  operator,  the  method  described  in  Lions  [11]  leads  to  an  integro- 
differential  equation  involving  e  . 

However,  one  of  the  principal  objectives  of  the  present  work  is  to  study 

properties  of  Galerkin  approximations  of  (3.1)  and,  in  particular,  to  obtain  a 
priori  estimates  for  such  approximations.  The  general  method  employed  by  Lions 

does  not  lead  to  results  from  which  a  constructive  approximation  theory  can  be 
easily  established.  For  this  reason,  we  give  here  an  alternative  proof  of 
Theorem  3.1.  which  was  suggested  by  Lions  [11],  also  based  on  the  notion  of 
elliptic  regularization,  but  in  which  a  different  form  of  the  regularized  problem 
is  used.  We  will  show  in  Section  9  that  the  constructive  nature  of  our  proof 
is  useful  in  studies  of  Galerkin  approximations.  Our  method  generalizes  those 
used  by  Lions  in  his  study  of  nonlinear  parabolic  problems  [10]  and  by  Dubinskii 
in  the  analysis  of  parabolic  problems  with  semibounded  variation  [5]. 

4.  An  Elliptic  Regularization 

We  begin  by  recalling  that  U  denotes  the  separable,  reflexive  Banach 
space  defined  by  (2.5)  which  is  everywhere  dense  in  III  .  We  will  denote  by 
[»,»]y  duality  pairing  on  U'  x  U  . 

We  next  introduce  a  family  of  "elliptic"  operators  A^:  U  -+  U' ,  e  is  a 


positive  real  number,  defined  by 
[Ac(u),v]u  - 


_9u  9v 
3t*3t 


H 


dv) 

U’^H 


+  [A(u),v] 


+  (u(T) ,v(T) ) 

u, v  e  u 


(4.1) 


where  A  is  the  operator  in  (3.1)  which  satisfies  conditions  (AI)  and  (All). 


The  problem  of  finding  u£  £  U  such  that 

[Ae(u£)  =  [f,v]  +  (uq,v(0))  V  v6  U  (4.2) 

is  an  elliptic  regularization  of  (3.1)  obtained  (formally)  by  adding  to 
3u/3t  +  A(u)  the  term  -e3^u/3t^  .  We  will  first  show  that  (4.2)  is  solvable 
and  then  prove  that  solutions  to  (3.1)  are  obtained  as  e  -*•  0+  . 

Lemma  4.1  The  operator  A£:  U  -*■  U’  defined  by  (4.1)  is  (i)  LI -pseudomono tone 
and  (ii)  coercive. 

Proof.  (i)  We  observe  that  the  operator  A£:  U  -*■  U'  can  be  expressed  by 
the  sum  A£  =  B£  +  A  where  B£  is  the  linear  operator  on  U  defined  by 


3u  3v1 

3v 

[3t’3t! 

H 

u*3Ej 

+  (u(T) ,v(T) )  ;  u,v  6  U  (4.3) 

We  prove  first  that  B£:  U  -*•  U’  is  positive  (monotone)  and  continuous.  Indeed, 
[Beu,u]u  -  e|u  jfj  +j  |u(0)  |  2  +  ±  | u(T)  |  2  >  0 

and 


[B£u,v] u  =  (u,v)H  +  (u,v)^  +  (u(0),v(0)) 

i  e|ulH  MW  +  lulH  |v|w  +  |u(0) |  | v (0) | 

1  (£  +  k1  +  k^)  HI  u i|| a  HI v||| u 

where  k^  and  k2  denote  the  continuous  embedding  constants  of  and 

U C( [0,T] ;H)  ,  respectively. 

Next  note  that  because  of  condition  (AI),  A  is  U-pseudomonotone  as  an 
operator  from  U  into  U'  .  Hence  the  operator  A£  is  pseudomonotone  since 
it  is  the  sum  of  a  continuous  monotone  linear  operator  and  a  pseudomonotone 
operator  (cf.  [11,  p.  189]). 


(ii) 

(A£(v)  .v^ 


I  v  (T) 


+  [A(v),v] 


4.3 

Hence 

Since  A  is  coercive  on  \)  by  assumption  (All),  the  first  term  in  the  last 
result  •+■  +  “  as  1 1 1  v | j|  y  -*■  00  and  this  proves  the  coercivity  of  A£  on  U  .  □ 
Therefore,  since  Lemma  4.1  establishes  sufficient  conditions  for  the 
surjectivity  of  A£:  U  ->  U'  (cf.  [11]),  the  following  existence  theorem  holds. 

Theorem  4.1.  For  any  fixed  e  >  0  ,  there  exists  at  least  one  solution 
u£  &  U  of  problem  (4.2). □ 

5.  Proof  of  Existence  Theorem  3.1 

We  now  return  to  the  existence  theorem  (Theorem  3.1)  for  the  nonmonotone 
parabolic  equation  (3.1)  in  which  the  operator  A:  [/-*■(/'  satisfies  conditions 
(AI)  and  (All).  Up  to  this  point,  we  have  shown  that  the  elliptic  regularization 
(4.2)  has  a  solution  in  U  for  any  £  >  0  .  We  now  show  that  solutions  to  (3.1) 
are  obtained  as  £  -*■  0+  .  This  is  accomplished  by  following  the  standard 
procedure:  (a)  establishment  of  a  priori  bounds  using  the  boundedness  and 

coercivity  of  A  ;  (b)  passage  to  the  limit  as  £  ■*  0+  ;  (c)  use  of  pseudomono tonicity 
arguments . 

Lemma  5.1.  Let  u£  6 .  U  be  a  solution  of  (4.2).  Then,  for  every  e  >  0  , 

there  exist  positive  constants  C^,  C^t  C^>  and  ,  independent  of  £  ,  such  that 

i)  llluElll  fC.ji.e.,  Hi u  HI  <  C  and 
LW  ex 

“>  *  T>T„ic 2 

Ul)  | ue <0)  |  <  C3  and  |  u^d)  |  <  C4 

Proof.  The  first  inequality  in  (i)  and  inequalities  (ii)  and  (iii)  follow 
from  the  boundedness  and  coercivity  of  A  .  Indeed,  replacing  v  in  (4.2)  by 

u  ,  we  have 

£ 


[A£(v),v]u  [A(v),v]+e|v|^  1  | v(0)|2+i'v(T)  | 2 

III v III (|  IIMII  +  |v[„  IIM||U 


5.2 


+\  |uc(0)|2  +|  I  U£  (T)  |  2  +  [A(ue),uc] 


=  [  f ,  u  ]  +  (u  ,u  (0)) 

t.  °  c  2 

i  III  fill,  lll-JII  +  V  |uo|2  +  |uE(0) 

2b 


where  b  denotes  the  Young's  constant.  Setting  b  =  1  produces  the  inequality 


[A(u  ),u  ]  .  |u 

— 6 — -  <  Ilf  III,  +  ?  — 
III  UE  III  -  2  IlluJ 


l“Elllicl 


and  by  setting  b  =  Ji  the  bounds  (ii)  and  (iii)  in  (5.1)  follow  from 
£|iElj;  +  \  l“e(0)|2  (T)  1  2 


5  [||l*(ue>  111.  +  Hit  HU  l||uj|  +  l«. 


To  prove  the  remaining  bound  in  (i)  ,  we  identify  a  solution  of  (4.2) 

with  a  distribution  u  6  V’ (0,T):V)  in  the  sense  that  u  satisfies  the 


distributional  equation 

-  cuc  +  u  ••  f  -  A(u£) 

Then,  proceeding  as  in  [10],  it  follows  that  in  H 


(5.2) 


-  cu  (0)  +  u  (0)  =  u 

e  e  o 


u£(T)  =  0 


and,  by  integrating  (5.2), 

u  ft)  =  fTie(t~X)/C  [ f (t)  -  A(r,u  (T))]dT 
t  JQ  t- 


from  which 


ll|«T/EH  !  Ill*  -  A(ue)||l, 
L1  (-T , 0) 


<  |||f  -  A(u  )  HI* 


This,  since  A  is  bounded,  completes  the  proof.  0 


Lemma  5.2.  Let  {u^JC  U  be  a  sequence  of  solutions  to  problem  (4.2) 
obtained  as  r.  ■*  0+  .  Then  there  exists  a  subsequence,  also  denoted  {u^}  » 
and  functions  u  6  W  and  X  6.  V’  such  that,  as  E  -+  0+  , 


weakly  in  W  ;  in  fact, 
weakly  in  1/  and 


weakly  in  M' 


weakly  in  H 


iii)  A(u£) 

iv)  ue(0) 

“E«> 


X  weakly  in  l/' 

u(0)  weakly  in  H 

u(T)  weakly  in  H 


Proof.  The  convergence  results  (i) ,  (iii)  and  (iv)  follow  from  Lemma  5.1 
the  fact  that  the  Banach  spaces  W  ,  V  and  H  are  reflexive,  and  the  boundedness 
of  A  .  To  establish  (ii) ,  note  that  from  Lemma  5.1  /e|u|  <  C„  ,  £>0. 

£  H 

Hence,  V  w  EH, 

•  j — ^  1  I 

£  (u  ,w)  <  / £  C  j  w  |  ->  0  as  £  *•  0 

e  U  1  H  □ 

Therefore,  by  virtue  of  Lemma  5.2  and  Green's  formula  (2.9),  we  see  that 
in  the  limit  as  £  ■*  0+  equation  (4.2)  becomes 

[u,v]  +  [ X » v ]  =  [f,v]  +  (uq-u(0) ,v(0))  V  v6  U 
Then,  using  the  fact  that  U  is  dense  in  U  and  property  (2.11),  we  conclude 
that  the  limit  u  of  u^_  satisfies  the  following  equation: 


|pv  +  [x,v]  =  [f,v]  V  v  6  1/ 


u(0)  =  u  in  H  J 

o 

It  remains  to  be  shown  that  x  =  A(u)  in  V'  •  Toward  this  end,  we  observe 
that  from  (4.2),  Lemma  5.2  and  (5.4) 

lim  { £  j  u  1 3  ~  [u  ,11  ]  +  (u  (T),u  (T))  +  [A(u  ),u  ]) 

ln  t 'n  l  c  e  c  c  c 


=  [f,u]  +  (uQ,u(0)) 

=  lim  {- [ u, u  ]  +  (u  (T),u(T))  +  [A(u  ) ,u ] } 
£  +  0 


5.4 


from  which  it  follows  that 


lim  [A(u  ),u  -u]  =  -lim  <  C | u  | ^ 

e+0  e+0  [ 


+  |  j u£ (0)-u(0)  |  2  +  \  |u£(T)-u(T)|2| 

1  0 

Hence,  since  A:  V  -*■  V’  is  W-pseudomonotone, 

lim  inf  [A(u  )  ,u  -v]  2.  tMu)  ,u-v]  V  v  6.  1/ 
eiO 

and 

lim  inf  [A(u  ),u  -v]  <  lim  sup  [A(u  ),u  ]  -  [x,v] 
eiO  e+0 

_<  lim  sup  [A(u  ),u]  -  [X.v] 
e+0 

=  [x.u-v] 

Consequently, 

[A(u)-x,u-v]  _<  0  V  v  6  V 
and  we  conclude  that 

X  =  A(u)  in  V’  (5.5) 

The  proof  of  Theorem  3.1  now  follows  immediately  from  (5.4)  and  (5.5). 

6.  Some  Sufficient  Conditions  for  Pseudomonotonicity 
We  will  review  briefly  here  some  results  of  Oden  [13]  which  provide  useful 
tests  for  pseudomonotonicity  of  a  certain  class  of  operators.  We  first  state 
a  corollary  of  Aubin's  compactness  theorem  [2]. 

Theorem  6.1.  Let  V  be  a  Banach  space  in  which  V  is  continuously 
embedded  and  consider  the  Banach  space  V  defined  by 


y  =  jv:vev,v  =  !^e  lPi(o,t;v1), 

1  <  Pi  1 00 1 


V  +  V 


>  (6.1) 


.  v  e  y 

'  Fi 

L  (0,T;V1) 

If  X  is  a  Banach  space  continuously  embedded  in  and  in  which  V  is 


mmKA 


«*vm- 


6.2 


compactly  embedded: 
c 

VQXQV1  (6.2) 

then  the  injection  from  V  into  the  space  LP(0,T;X)  is  compact: 
c 

VQLP(0,T;X)  <6‘3) 

□ 


Following  Oden  [13],  let  A(u)  denote  values  of  an  operator  A  from  the 
Banach  space  V  into  its  dual  V'  which  has  the  property  that  there  exists  a 
map  (u,v)  H-  A(u, v)  ,  t  6  (0,T)  ,  from  V  x  V  into  l/’  such  that  A(u,u)  =  A(u) 
and  the  following  conditions  hold: 

BI.  V  v  6  V  ,  u  A(u,v)  is  hemicontinuous  from  1/  into  l/'  ;  i.e., 

V  u,v,w  6  1/  ,  the  function 

<J)(s)  =  [A(u+sw,v)  ,w]  ,  s  6,  R 
is  continuous  in  s  . 

BII .  V  u,v€B^(0)  »  {w  6  V:  IIHH  <  W,  U  >  0}  , 

[A(u,u)-A(v,u)  ,u-v]  >  -H(y,  ||  u-v  ||  )  (6.4) 

1/(0, T;X) 

-f-  ^ 

where  H:  R  x  R  ■>  R  (R  =  [0,°°))  is  a  function  continuous  in  each  of  its 
arguments  with  the  property  that 

lira  4  H(x,0y)  =  0  V  x,y  £  R+  (6.5) 

0+0 


and  LP(0,T;X)  is  as  in  Theorem  6.1. 

Bill.  If  u  — *  u  weakly  in  V  of  (6.1),  then 


lim  inf  [A(v,u  )  -  A(v,u),u  -u]  _>  0  V  v  £  1/ 

n-*» 

and 


\  (6.6) 


lim  inf  [A(v,u  )  -  A(v,u),w]  =  0  V  v,w  6.  V 

n-xio 

BIV.  A:  V  ■*  V1  is  bounded. 

The  importance  of  these  conditions  is  made  clear  in  the  following  theorem 
proved  in  [13]. 

Theorem  6.2.  Let  A:  V  ->■  V'  satisfy  conditions  (BI) ,  (BII),  (Bill),  and  (BIV). 


Then  A:  V  ■*  V’  is  V-pseudomonotone.  0 
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PART  II.  A  MODEL  PSEUDOMONOTONE  DIFFUSION  PROBLEM 
7.  Existence  Analysis 

In  this  part,  we  are  concerned  with  the  qualitative  analysis,  Galerkin  and 
Faedo-Galerkin  approximations  of  the  following  generalized  nonlinear  diffusion 
problem:  Let  ft  be  an  arbitrary  bounded  domain  in  Rn  with  boundary  9ft  , 
and  0  <  T  <  00 

find  u  =  u(x,t)  ,  (x, t)  £  ft  *  (0,T)  ,  such  that 


Given  data  f  in  ft  *  (0,T)  and  initial  data  uq  on  ft  , 


-  V  •  a(Vu)  +  b (u,Vu)  =  f  ,  in  Q  =  ft  x  (0,T) 

u=0,  on  Z=9ftx  (0,T) 


u(*  ,0)  =  u  ,  on  ft 


(7.1) 


where 


<l(Vu)  =  aVu  +  k|Vu 


P-2 


Vu 


a,k  €  L  (ft)  ,  2  <  p  < 


a(x)  _>  aQ  _>  0  and  k(x)  _>  kQ  >  0  a.e.  x  £  ft 


(7.2) 


and  b(u,7u)  =  b(x,u(x,t),  Vu(x,t))  is  subject  at  a.e.  (x,t)  6  Q  to  the 
conditions : 


Cl.  |b  (Ctp 

1  <c|d'  |j|c  ,  V 

(C.O  6  R  : 

<  Rn 

q  =  0  or 

q  >  1  ,  r  =  0  or 

r  2 

1 

\ 

> 

(7.3) 

l_<q+r<p-l 

( 

i 

CII .  b  (?,p 

is  totally  Frechet  differentiable  in  IR  x  (Rn 

and 

its  partial 

derivatives  9^b: 

R  x  Rn  -*•  L  (R,R)  and 

:  R  x 

Rn  -*■  L  (Rn,R) 

are 

such  that,  for 

(q,r)  satisfying 

(7.3)  and  V  (£,£)  6  R  * 

Rn  , 

|3cb(C.C) 

1 1  c,  ur1  ujr . 

if 

q  i  0 

|9rb(C,0 

i  <  cr  id’  ur1 . 

if 

r  4  0 

The  case  r  =  0  will  be  understood  as 

b  = 

b  (u) 

(not  function  of 

Vu) ,  and  the 

case  q  =  0  as  b  =  b(Vu)  (not  function  of  u) . 

In  this  case,  we  take  as  spaces  V  and  H  the  usual  Sobolev  spaces 


7.2 


V  =  W1,P(52)  ,  2  <  p  <  00 

O  — 

H  =  L2(52) 

Then,  with  the  conventions  of  Section  2  in  force,  the  model  problem  (7.1)  assumes 
the  following  form: 

Find  u  6,  W  such  that 


(7.4) 


+  A(u)  =  f  ,  f  given  in  U' 

2 

u(0)  =  u  ,  uq  given  in  L  (52) 
where  A:  V  ■*  V'  is  defined  by 

[A(u),v]  =  [A1(u),v]  +  [A2(u),v] 


(7.5) 


(A1(u) ,v] 
[A2(u) , v ] 
in  which  a(\7u) 


=  a(x,Vu(x, t)) *Vv(x,t)  dxdt 

Jq  ~ 

* 

=  b (x,u(x, t) ,Vu (x, t))v(x, t)  dxdt 

JQ 

is  as  defined  in  (7.2)  and  b(u,Vv) 


(7.6) 


is  subject  to  conditions 


(Cl)  and  (CII). 


We  now  proceed  to  establish  the  existence  of  solutions  to  problem  (7.5). 

The  following  two  theorems  determine  fundamental  properties  of  the  operator  A  . 

Theorem  7.1.  Let  A:  1/  ■*  V1  be  the  operator  defined  in  (7.6).  Then 
i)  A  is  bounded,  ii)  A  is  coercive,  and  iii)  A  is  locally  Lipschitz 
continuous  in  the  sense  that  V  u,v  €.  Bp(0)  =  {v  6  1^:  Hi v III  <  y  ,  y  >  0}  , 
w  &  1/  ,  there  is  a  positive  constant  C(y)  such  that 

|[A(u)-A(v),w]|  <  C(vt)  |||u-v|||  lj|w|||  (7.7) 

Proof.  We  shall  use  the  notation  a  =  j|a||  , 

l”(52) 


k  = 


L  (52) 


and 


II- 


II  s 

L  (Q) 


i)  Applying  Holder's  inequality,  we  obtain,  V  v,w  6  V 


I A^  (v)  ,w]  |  jc 


'  aJVvl  +  kJVv!P~1J  |Vw|dQ 

ammes(Q)(P  2)/P||  Vv||  p>Q  +  k^ 


P-1 

P.Q 


Vw 


P.Q 


7.3 


|  £  A  (v)  ,w]  |  _<  c  |  v ) q  |  Vv  1 r  j  w  |  dQ 


<  c  mes  (Q^P_1  q"r)/p||  v||  q>Q||  Vv||  p>Q||  w||  p>Q 
Hence,  p-2 

iii*(»)in,  p  111 v  111  ♦kjMir1 

p-l-q-r 

+  c  mes (Q)  P  |||vj||q+r  ,  V  v  6  V  (7.8) 

ii)  From  Friedrichs'  inequality,  it  follows  that,  V  v  g.  LS(0,T;W^,S(fl))  , 

1  <  s  <  00  , 

II  vll  Ss  !  s  *  II  v  II  s  Q  +  II  Vv!l  s  Q 

LS(Q,T;Wi,S(Q))  S>Q  S‘Q 

o 


<  (cS(s,n)  mes  (ft) S /n  +  1)  ||  Vv||  ®  (7.9) 

—  S 


Thus , 


[A(v),v]  > 


(a  |Vv|2  +  k  |Vv|P)dQ  -  c | v | q+^ | Vv| rdQ 

Q  °  °  Q  _! 

2  aQli  VvJ|  2^  +  kQ(l  +  cP(p,n)  mes(ft)p/n)  |||v|||P 
p-l-q-r 

-  c  mes (Q)  P  |||v|||1+q+r  ,  V  v  6.  V  (7.10) 

But  using  Young's  inequality  (2.12)  in  the  last  term  in  (7.10)  and  choosing  b 
small  enough  leads  to 

[A(v),v]  >  aQ]|  Vvj|  2>Q  +  T1ll|v|HP  -  Y2T  ,  V  v  6  1/  (7.11) 

where  Y^  and  Y2  are  >  0  .  Therefore, 


lA(v.)  >.y!  >  „  111..111P-1  t2T 


Y  v 


•>  +  00  as  v  00 


iii)  By  the  inequality  in  R  [15] 

llx|r~2  x  -  |y|r  2y|  1  c{|xj  +  |  y  [ }  r~ 2 1  x-y  | 
c  =  /r  -  T  if  2<^r_<3,  c  =*  r  -  1  if  3<r<°° 
and  Holder’s  inequality,  we  obtain,  V  u,v  6  B  (0)CZ  1/  ,  w  61/  , 

| [A1(u)-A1(v) , w ] | 

[aJV(u-v)|  +  kmc(p)(|Vu|  +  |  Vv  j  )P_2|V(u-v)  |1  |Vw|dQ 

t  B=A  |  J 

mes(Q)  p  +  k  c (p)  (2p)P_2  >  ||| u-v|||  |||w||| 


(7.12) 


(7.13) 


1 

i 


7.4 


We  now  use  hypothesis  (CII).  First  we  observe  that 

U 


lIJfpUUUCOAO  \  »» w  wwww*.  *  V.  % 

[A2(u)-A2(v),w]  =  |  j*  -db^7-^-  w  d6dQ 


{3^b(C,VOn  +  3v^b(C,V0-Vn}w  d0dQ  (7.14) 

where  £  =  v  +  0r|  ,  rj  =  u-v  and  0  6  [0,1]  .  Hence,  because  of  (CII)  and  Holder's 
inequality, 

| [A2(u)  -  A2(v)  ,w  ] | 

1  ff  (c  Ulq_1|vc|r|n|  +  c  Ulqlvc|r~1|vn|]  |w|dQd0 

JqJq  t  q  J 


<  (c  +  c  )  mes  (Q) 
—  q  r 

Then,  since  u,v(  B  (0)  Cl  V  , 


EiSZE  ! 

"  ! 
'0 


l5lll,+r"1de]||n| 


P.Q 


[A2(u)-A2(v),w]  |  <  y3yq+r  1|||u-v[ 


w 


LP(Q) 


p-l-g-r 


y  =  (c  +  c  )  mes (Q) 
J  H 


(7.15) 


Therefore,  from  estimates  (7.13)  and  (7.15),  (7.7)  follows  and  this  completes 
the  proof  of  the  theorem.  Q 

The  next  property  of  A  ,  established  below,  is  crucial,  not  only  in  proving 
the  existence  of  solutions  to  (7.5)  but  in  subsequent  studies  of  approximations. 
Theorem  7.2.  The  operator  A:  V  V’  defined  in  (7.6)  satisfies  the 

O 

following  nonlinear  Garding-type  inequality: 


[A(u)-A(v),u-v]  >  a  a  ||  u-v||  2  .  +  a,  |||u-v||| 

0  °  L  (0,T,H^ (H) ) 

-  a2(y)||  u-v ||  p 

L  LP(Q) 

V  u,v&B^(0)  =  (w6  l/:  j||w|||  <  y  ,  y  >  0} 

1  12 

Here  Hq((5)  =  Wq'  (fi)  and  aQ  ,  ,  a2(y)  are  constants  satisfying 

aQ  >  0  ,  a2  >  0  ,  a2(y)  >  0  . 

Proof.  We  observe  that,  V  u,v  6  V  , 


(7.16) 
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( A(u)  -  A(v),  u  -  v]  j>  [A1(u)  -  A1(v),  u  -  v] 

-  | [A2(u)  -  A2(v) ,  u  -  v] | 

From  the  inequality  in  [15] 

(|x|r  2x  -  | y  | r  2y,x-y)  _>  21  r|x-y|r  ,  2  _<  r  <  ° 

and  (13.9),  it  follows  that 
[A^^  (u)-A1  (v),u-v] 

>  jaj  V(u-v)  |  2  +  ko21_P|'7(u-v)  |PJ  dQ 


(7.17) 


(7.18) 


>  a 
—  o 


l+c2(2,n)  mes(ft)2/n 


-1 


u-v 


L2(0,T;W*,2(ft)) 


+  k  2 
o 


1-p 


l+cP(p,n)  mes(fi)p/,n] 


-1 


(7.19) 


On  the  other  hand,  according  to  (7.15)  and  Young's  inequality  (2.12),  for 

u,v  6  B  (0)C  V  , 

^  P 

|[A2(u)-A2(v),u-v]|  £  ~  |||u-v|||P 

p'(q+r-l)p’ 


y3  U 


p’bP’ 


u-v 


LP(Q) 


(7.20) 


Therefore,  introducing  (7.19)  and  (7.20)  into  (7.18)  and  choosing  b  small 
enough,  the  desired  result  (7.16)  is  obtained.  D 

2 

Theorem  7.3.  For  any  data  f  6  V'  and  uq  6  L  (ft)  ,  there  exists  at  least 

one  solution  u  to  problem  (7.5). 

Proof.  Theorem  7.1  confirms  that  conditions  (BI),  (Bill)  with  (/  =  W  ,  and 

(BIV)  of  Theorem  6.2  are  satisfied  ((Bill)  is  trivially  satisfied).  Condition  (BII) 

with  X  =  LP (ft)  (QV1  and  in  which  V  is  compact  (cf.,  e.g.,  [1]))  also  holds 

since,  by  virture  of  (7.16), 

[A(u)-A(v)  ,u-v]  _>  — H (n ,  | j  u-v||  ) 

LP(Q) 

V  u,v  B^(0)C1  V  ,  where 

H(y.y)  =  a2(y)yp  ,  y£R+  ,  p’  =  p/p(p-l) 

Therefore,  A  is  coercive  and  W-pseudomonotone  from  V  -*■  l/1  and,  by  virtue  of 


7.6 


Theorem  3.1,  the  assertion  of  the  theorem  follows.  0 

Remark  7.1.  From  the  proofs  of  Theorems  7.1  and  7.2,  it  is  apparent 

that  the  operator  A  of  (7.6)  regarded  as  a  map  from  V  into  V'  ,  is  bounded, 

© 

coercive  and  locally  Lipschitz  continuous,  and  satisfies  the  Garding-type  inequality 
<£(u)-A(v),u-v>  >  a  a  ||  u-v||  2  +  a  ||  u-v||  P 

-  a,(p)||  u-v||  P’  ,  V  u,v  6  B  (O)CV 

J  LP(ft)  P  (7.21) 

According  to  Oden  [13],  A:  V  ■+  V'  is  necessarily  (V-) pseudomonotone.  Hence, 

from  the  theory  of  pseudomonotone  elliptic  equations  (cf.  [11]),  A  is  surjective 

from  V  -+■  V’  ;  i.e.,  there  exists  at  least  one  solution  in  V  to  the  stationary 

problem 

A(u)  =  f  ,  f  given  in  V1  (7.22) 

The  evolution  problem  (7.5)  possesses  at  least  one  equilibrium  state  for  each 
f  £V  .  □ 


Sufficient  Conditions  for  Uniqueness 


We  now  proceed  to  determine  sufficient  conditions  for  uniqueness,  of 


solutions  to  the  pseudomonotone  diffusion  problem  (715). 

In  the  case  of  monotone  parabolic  problems,  "monotonicity"  "uniqueness" 

and  this  follows  from  the  differential  inequality  of  Caratheodory  type: 

d|u(t)-v(t)|  /dt  _<  0  ,  a.e.  t  £ [0,T]  ,  |u(0)-v(0)|  =  0  ,  whose  unique  solution 

2 

is  |u(t)-v(t)|  =  0  ;  u  and  v  are  supposed  solutions  of  the  problem.  This 

O 

suggests  that  in  the  non-monotone  case  with  Garding-type  inequality,  the  possibility 
of  establishing  a  differential  inequality  of  the  form 

~  |u(t)-v(t) | 8  <  a|u(t)-v(t)|s 

a6R,  2  <  s  <  ®  ,  for  a.e.  t  £  [0,  T]  *  (8.1) 

|u(0)-v(0) | 8  =  0 


8.2 


would  be  sufficient  for  concluding  uniqueness.  Indeed,  from  Theorem  3  of  Olech 
and  Opial  [14],  |u(t)-v(t)|S  =0  is  the  unique  solution  to  (8.1).  We  show  that 

in  certain  particular  cases  and,  in  general,  for  sufficiently  smooth  solutions 
of  problem  (7.5),  this  is  the  case. 

Theorem  8.1.  Let  u  6  W  be  a  solution  of  problem  (7.5).  Then  u  is 
unique  in  the  following  three  cases: 


i)  r 

=  0 

and 

q  =  1 

ii)  r 

-  0 

and 

n  <  p 

> 

(8.2) 

iii)  a 

o 

>  0 

and 

u  6  L  (0,T;W^  (O')) 

Proof.  Assume  that  u  =  u(t;f,uo)  and  v  =  v(t;f,uQ)  are  two  solutions 
of  problem  (7.5)  and  define  r)  =  u-v  .  As  is  apparent  from  (7.19), 

o''  _ 

^  ,  ..  II  _  II  p 


<^A1(u(t))-A1(v(t)),n(t5>  >  aoao||  n(t)  ||  ±  +  ajl  n(t) 

H  (H) 
o 

for  a. e.  t  6  [0,T] 


(8.3) 


where  aQ  >  0  and  >  0  .  Thus,  from  the  difference  of  the  equations  satisfied 
by  u  and  v  ,  we  obtain  the  integral  inequality 

'  fT 


I  In(x)|2  +  aoao 


0 

rT 


n(T)||  .  dt 
IT(fi) 


<k  (u(t))-A  (v(t)),n(t£>  dt 


a  >  0  ,  V  re  [0,T] 


(8.4) 

We  now  estimate  the  right-hand  side  term  via  the  formula  (7.14)  with 


w  =  n 


i)  r  =  0  and  q  =  1  .  In  this  case  we  have  the  estimate 


<X  (u(t))-A  (v(t)),n(t£>  dt 


<  c 


|n(0|  dt  , 


vie  [o,t] 

which  combined  with  (8.4)  gives  the  integral  inequality 


(8.5) 


*In  this  case,  the  question  of  existence  appears  to  be  open. 


8.3 


|n(T)|2  <  2c 


fT 


I n (t) |  dt  ,  Vie  [0,T] 


(8.6) 


But  this  is  equivalent  to  (8.1)  with  a  =  2cq  and  s  =  2  .  Consequently,  n  =  0  . 

ii)  r  =  0  and  n  <  p  .  From  the  Sobolev  embedding  theorem  (cf.  [1,  Chap.  5]), 
W*,P(ft)  is  continuously  embedded  in  C_(ft)  =  {v  6C (ft):  v  bounded  in  0}  whenever 

O  15 

n  <  p  .  Then 

V  =  LP(0,T;wl’P(ft))QLP(0,T;L°°(ft))  ,  n<p  (8.7) 

Let  y  be  chosen  such  that  u,v  6  B^(0)C1  1/  .  Then,  from  (7.14)  with  w  =  n 
and  using  (8.7),  we  obtain,  V  x  6  [0,T]  , 

ft 

<A  (u(t))-A  (v(t)),n(t£>  dt 


2  /  “2 

0  flfT 

< 


0J0J 

rlfx 


_  1  ry 

c^|g(x,t)|  |n(x,t) I  dxdtdQ 


T 

ojo 


ca  II  5(011  ^  |n(t)  I 2dtd0 

q  L  (ft)  2 

ft 


<  c  y 

-  q 


q-l 


I n(t) |  dt 


2  <  s 


=  2P 
p+l-q 


<  P 


(8.8) 


Introducing  this  estimate  into  (8.4)  produces  the  integral  inequality 
(  1'1s/2fT 

,s  -  L  q_1l  I  I n (t) I sdt  ,  V  x  e  [0,T] 


I n (x)  |  < 


2c  y 

q 


(8.9) 
q-l 


2c  y 

q 


which  is  equivalent  to  the  differential  inequality  (8.1)  with  a  ■ 

and  2  <■  s  *  sp/ (p+l-q)  <  p  .  Therefore,  y  =  0  . 

iii)  aQ  >  0  and  u  6  L  (0,T;W^’  (ft))  .  Let  y  >  0  be  such  that 

u,v  6  B^(0)C1  L  (0,T;W^*  (ft))  .  Then,  from  (7.14)  with  w  =  n  ,  we  obtain 

V  x  6  [0,T] . 

rx 

<A2(u(t))-A  (v(t))  ,n(t)^>  dt 


|v^(x,t) |r|n(x,t) | 


+  cr|C(x,t)|4|VC(x,t)|r"1|Vn(x,t) 


q+r-l 


^c^c(2,n)  mes(ft)^n+cr 


n(t) 


|r|(x,t)  |dxdtd8 
I n ( t) I dt 


H  (ft) 
o 


(8.10) 


i/2 


ispop  -muuuw 


9.1 


Hence,  since  by  hypothesis  aQ  >  0  ,  we  can  apply  Young's  inequality  with 
constant,  e.g.,  b  =  ,  to  obtain  the  upper  bound  for  (8.10) 

a  a  fT  „  .  cT 


— [  ||  n(t)  ||  2  dt  +  f(  |n(t)  |2dt 
J0  n  (Q)  ZJn 


where  a  =»  a(l/b  )  >  0  .  Combining  these  results  with  (8.4)  gives 


9  fT 

|n(x)  |  <_  a 

Jo 


n(t) |  dt 


(8.11) 


and,  consequently,  (8.1)  holds  with  s  =  2  and  n  =  0  .  □ 

9.  Galerkin  Approximations 

In  this  section,  we  study  Galerkin  approximations  of  the  model  problem 
(7.5)  which  are  based  on  an  elliptic  regularization  of  (7.5)  obtained  using 
the  ideas  described  in  Section  4.  We  will  establish  some  results  on  the  strong 
convergence  of  such  approximations. 

For  the  model  problem  (7.5),  we  introduce  the  corresponding  elliptic 
regularization  (4.2): 


e(u£»v)^  -  (u,,v)H  +  (u£(T),v(T))  +  [A(u£),v] 
=  [f,v]  +  (uq,v(0) )  ,  ¥  v  £  U 


(9.1) 


where  A  is  the  operator  defined  in  (7.6).  According  to  Theorems  7.1  and  4.1, 
a  solution  u£  6  U  exists  to  such  a  regularization  for  every  e  >  0  .  Moreover, 
according  to  Section  5  ,  in  the  sense  of  V'  C.  L(P((0,T)),W  ((2))  ,  u£  satisfies 

the  distributional  equation 


-eii  +  u  +■  A(u  )  =  f  in  V1 

E  E  £ 

-eu£(0)  +  u£(0)  =  uq  in  L2(£2) 
u£(T)  =  0  in  L2(fi) 


(9.2) 


which  is  equivalent  to  (5.1),  and  for  any  sequence  {u  }  C  U  of  solutions,  there 

e  £>0 

exists  a  subsequence,  also  denoted  ^U£^£>o  such  that,  as  e  -*■  0+  ,  u£  converges 
weakly  to  a  solution  u  of  (7.5)  in  the  sense  of  Lemma  5.2  with  X  “  A(u)  . 


To  construct  Galerkin  approximations  of  (9.1),  we  introduce  a  family  of 


9.2 


subspaces  {U  }  of  U  such  that  i)  U.  is  finite-dimensional  with 

h  0<h<l  li 


basis  functions  {0^,  0^* 


“h 


}  ,  with  dimension 


-*■  00  as  h  -*•  0  and 


ii)  U  U  is  dense  in  U  .  A  Galerkin  approximation  of  (9.1)  involves  seeking 
h  h 

a  function  U  €  U,  such  that 
£  h 


+  t»(^>,*k] 

=  [f,0kl  +  (uo.4>k(0))  ,  k  =  1,2,...  ,m^ 


(9.3) 


The  solvability  in  of  (9.3)  is  assured  by  Lemma  4.1.  Similarly  as  in 

the  proof  of  Lemma  5.1,  if  {uj?}(}<h<l  is  a  se9uence  of  Galerkin  approximate 
solutions,  it  can  be  shown  that  there  exist  constants  K^,  K2,  and  K^,  independent 

of  h  ,  such  that  ll|U^  III  i  K1  •  l^l«  i  K2  •  lU£<°>  l  5  K3  and  |^(T)|  £  K4  . 

Then,  via  weak  compactness  and  pseudomonotonicity  arguments  (as  those  used  in 
Section  5),  it  follows  that  there  exists  a  function  u^.  and  a  subsequence,  also 
denoted  »  such  that,  as  h  0+  , 


^“"U£ 

weakly  in 

V 

•  h  _ v  • 

U£-^U£ 

weakly  in 

l2(q) 

A(U^)— *  A(u£) 

weakly  in 

V 

U^(O) — ^ue(0) 

weakly  in 

L2(ft) 

U^(T)— ^ue(T) 

weakly  in 

L2(ft) 

(9.4) 


We  will  now  demonstrate  that  for  our  model  problem  (7.5)  much  stronger 
results  can  be  obtained. 


Theorem  9.1.  Let  U  be  a  weakly  convergent  subsequence  of 

solutions  to  problem  (9.1)  and  let  u&  W  be  the  corresponding  weak  limit,  solution 
of  problem  (7.5).  Then,  as  £  -*■  0+  , 

strongly  in  1/ 


u  -*  u 

£ 


71  u 


0 


u  (0)  -*■  u 
e  o 

ue(T)  -  u(T) 


strongly  in  L  (Q) 
2 

strongly  in  L  (12) 

2 

strongly  in  L  (Q) 


(9.5) 


9.3 


Proof.  We  regard  equation  (7.5)  as  holding  on  U  and  suits  t  met  (‘>.1) 

from  it.  The  following  orthogonality  condition  is  obtained: 

-(Cu  ,v)H  +  (u  -u  (0) ,v(0)) 

£  rf  o  £ 


+  [u-u^,v]  +  [A(u)-A(u^) , v]  =0  V  v  €  U 


(9.6) 


According  to  the  a  priori  bound  (5.1)^,  there  is  a  y  >  0  independent  of  £  , 

O 

such  that  u^,  u  €.B^(0)C  ^  •  Hence,  using  formula  (2.9)  and  the  Garding-type 

inequality  of  (7.16),  we  see  that 

Iuo~ue(0)|2+  [u— u£,u-u£]  +  [A(u)-A(uc) ,u-u£] 

>\  i u0_u£ (0) i 2  +  ~  |u(T)-uc(T)|2 

+  “illlu-^lll15  “  a2(p)||  u-u  ||  P 

1  £  Z  LP(Q)  (9.7) 


Next,  combining  these  two  results,  we  conclude  that 

\  |uo-u£(0)|2  |u(T)-ue(T)|2  +  aJIu-uJH15 

<_  a  (y)  ||  u-u  ||  P  +  (u  — u  (0)  ,u  — v(0) ) 

1  LP(Q)  °  C  ° 

+  [u-uc  ,u-v]  +  [A(u)-A(uJ  ,u-v]  +  (eu^.v)^  -  1 71  u 


V  v  6  U 


(9.8) 


Due  to  the  compact  embedding  of  W  in  LP(Q)  (cf.  Theorem  6.1)  and  the  weak 

convergence  result  of  Section  5,  (9.5)  follows.  □ 

Theorem  9.2.  Let  (Uh  €.  U,  ,  be  a  subsequence  of  Galerkin  approximate 

-  £  h  U<h<i 

solutions  defined  by  (9.3),  converging  weakly,  in  the  sense  of  (9.4),  to  a  solution 


u£  &  U  of  problem  (9.1).  Then,  for  fixed  £  >  0  ,  as  h  -*■  0 


l/1  ■*  u  strongly  in  V 


U*1  ->  u 

£  £ 


strongly  in  L  (Q) 


U^(0)  •*  u^(0)  strongly  in  L2(fi) 
U^(T)  -*■  u^(T)  strongly  in  L2(f2) 


Proof.  We  follow  similar  arguments  to  those  given  previously.  Restricting 
(9.1)  to  and  subtracting  (9.3)  from  it,  we  obtain  the  orthogonality  condition 


9.4 


+  [A(u  )-A(U'  ),W]  =  0  ,  V  W  6U, 

e  e  n 


e(u£-uj!,W)H  -  (u£-U£,W)H  +  (u£(T)-l£(T),W(T)) 

(9.10) 

Now,  from  (5.1^  and  since  III III  Ki  v  h  €  (0,1]  there  is  a  y  >  0  , 
independent  of  h  such  that  U£,u£  €  B^(0)CL  V  .  Then,  by  virture  of  (2.9)  and 
(7.16),  it  follows  that 

Eiv“? « -  +  K<™£<t)I2 

+  [A(u£)-A(U£) ,U£-U£] 

>e|;e-6^+i  |uc<0)-^(0)|2+{  |uE(T)-^(T)|2 

+  o  l|  u  -UE||  2  +  o  |||u  -UE|||P 

L  (0,T;H2(£2))  1  E  £ 


-  cc  (M)  II  u  -u»!|| 

C  £  LP(Q) 

Therefore,  combining  (9.10)  and  (9.11) 

£lv^e  !w  +  2  lue  (O)-U^(O)  |  2  +  \  |u.  (T)-U£ (T) | 2 


(9.11) 


+  a  u  -U 
o 


e-cH22  1  ♦“illlV^lil’’ 

LZ(0,T;H^(n)) 


i“2(p,IIV<ll^(Q)  +  E(V“e'V">« 

-  (ue-U£,u£-W)H  +  (u£(T)-U£(T),u£(T)-W(T)) 

f  [A(u  )-A(Uh),u  -W]  V  W  6U, 

c  c  e  h 


(9.12) 


But,  according  to  Theorem  6.1,  U  is  compactly  embedded  in  LP(Q)  and  U£ 
converges  weakly  to  u£  in  the  sense  of  (9.4).  Hence,  the  right  side  of  (9.12)  -*■  0 
as  h  *  0+  and  this  proves  the  theorem.  Q 

We  next  give  an  error  estimate  for  the  Galerkin  approximations  of  the 
regularized  elliptic  problem  (9.1). 


Theorem  9.3.  For  fixed  e  >  0  ,  let  u£  6  U  be  a  solution  of  problem  (9.1) 
which  is  the  strong  limit  (in  the  sense  of  (9.9))  of  a  subsequence  of  Galerkin 
approximate  solutions  {U£  ^  uh^0<h<l  defined  bF  (9.3).  Then  the  following 


9.5 


approximation  error  estimate  holds  V  W6U^  : 

\  C1|u.(0)-l£<0)|2  +  \  |uc(T)-U^(T)|2 

■T=Hi 

i  M  V"e"  "p(g)  +  c2l<*c(0)-w<0)|2  +  c3|ue-»|J 

+  C|||VU|||P’  +  C4|ic-Wj2  (9.13) 

where  C^,  i  =  l,...,4,ao>  6^  =  (^(a  ),  a2  =  a2(T,y)  e  =  e(e)  and  C  =  C(C(T,y)) 

are  strictly  positive  constants.  Here  C(T,y)  is  the  local  Lipschitz  continuity 
constant  of  (7.7). 

Proof .  The  estimate  (9.13)  follows  directly  from  (9.12)  upon  applying 
formula  (2.9),  the  local  Lipschitz  continuity  of  A  ,  (7.7),  and  inequality  (2.12).  □ 

10.  Faedo-Galerkin  Approximations 

We  are  concerned  here  with  Faedo-Galerkin  approximations  of  the  model 

pseudomonotone  diffusion  problem  (7.5).  We  note  that  this  type  of  approximation 

process  is  not  necessarily  well-defined  for  non-monotone  parabolic  problems:  the 

corresponding  weak  convergence  is  a  conditional  property.  We  shall  show  that 

Faedo-Galerkin  approximate  solutions  to  problem  (7.5)  exist  and  are  unique,  and 

determine  sufficient  conditions  for  weak  and  strong  convergence. 

Let  {V  }g<h<i  be  a  family  of  finite-dimensional  subspaces  approximating 

the  space  V(=  W^’^(f2))  in  the  following  sense:  (i)  }  denotes 

+  l  I  ^ 

a  basis  for  ,  with  dimension  m^  °°  as  h  -*■  0  ,  (ii)  S'  is  dense  in 

V  .  A  Faedo-Galerkin  approximation  in  of  problem  (7.5)  is  defined  as  an 

h  u 

absolutely  continuous  function  U  from  [0,T]  -+•  V,  ,  i.e.,  u  €  C.([0,T];V.)  , 

h  Ah 

solution  of  the  system 

<Uh(t),^k>+<A(Uh(t)),1pk>=<f(t),^k>  k  =  1,2 . mh  9| 


+  a 


(fl» 


|u  -u"| 
1  e  e 1 


+  c 


IvuJh 


+  a  a  u  -U  , 
o  o 11  t  c 


,.h  ii  2 


2 

L  (0 


Uh(0)  -  uh 

o 


(10.1) 


10.2 

h  2  i 

for  a.e.  t  6  [0,T]  and  where  Uq  -*■  uq  strongly  in  L  (Q)  as  h  -*■  0  .  We 

Jl  •  ll 

observe  that  if  U  is  solution  of  (10.1),  then  its  time  derivative  U  belongs 

to  1,P  (0,T;V  )  but  not  necessarily  to  l/'  =  LP  (0,!;W  (il))  (C(J  LP  (0,T;V,))  . 

n  ,  h 

h 

We  next  establish  the  solvability  of  problem  (10.1). 

Theorem  10.1.  For  each  h  6  (0,1)  ,  the  Faedo-Galerkin  approximation  problem 

(10.1)  possesses  a  unique  solution  U^*  £  C^([0,T];V^)  continuous  with  respect 

to  Uh  . 
o 

Proof.  The  local  existence  of  solutions  to  (10.1)  in  C.(f0,t,  1;V,  ), 

A  h  h 

t,  >  0  ,  is  implied  by  the  pseudomonotonicity  property  of  A  (cf.  Remark  7.1). 
h  ' 

Indeed,  f  £  V  and  A  is  necessarily  bounded  and  demicontinuous  from  V  -*■  V' 

and  these  are  sufficient  conditions  for  the  vector  field  F(t,U)  =  (<^f(t), 

m  m, 

<A(U(t)),i|ik>)  from  D  =  [0,T]  x  R  n  -  R  1  to  satisfy  the  Caratheodory 

“h 

conditions  in  D  .  Here  U  £.  IR  denotes  the  coordinate  vector  of  U6V, 

~  h 

with  respect  to  the  reciprocal  basis  of  V.  . 

h 

The  uniqueness  and  continuous  dependence  on  the  initial  data  of  local  solutions 
to  problem  (10.1)  follows  from  the  condition  [7]:  for  each  compact  set  w£D  , 
there  is  a  function  £.  L^(0,T)  such  that 

|F(t,U)-F(t,W)|  <  gw(t)|U-W|  ,  (t,U),(t,W)  e  w  (10.2) 

which  is  satisfied  because  A  is  locally  Lipschitz  continuous  from  V  -*■  V’ 

(cf.  Remark  7.1). 

It  remains  to  be  proved  that  the  interval  of  existence  [0,t^]  =  [0,T]  . 

This  is  a  consequence  of  the  coercivity  of  A  from  V  -*•  V'  ,  as  follows  from 
part  (1)  of  the  proof  of  Theorem  10.2  given  below.  □ 

We  now  proceed  to  analyze  the  convergence  of  the  Faedo-Galerkin  approximation 
process . 

Theorem  10.2.  From  the  sequence  of  Faedo-Galerkin  approximate  solutions 

defined  uniquely  by  (10.1),  there  is  a  subsequence,  also  denoted  {U^}rt  , 

0  <h<l 


and 
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there  exist  functions  u  €  U/  and  X  6  V'  such  that,  as  h  -*•  0  , 


weakly  in  V 


U  - Au  weakly*  in  L  (0,T;L  (ft) ) 

aoj11) - weakly  in  V 

U^T)  - u(T)  weakly  in  L2(ft) 


(10.3) 


3u 

1>’VJ 


+  [X.v]  =  [f.v]  ,  V  V  e  V 
u(0)  =  u 


(10.4) 


Moreover,  the  limit  function  u  is  a  solution  of  problem  (7.5)  (i.e.,  X  =  A(u)) 


provided  one  of  the  following  conditions  is  satisfied: 


i)  Uh  6  V'  ,  0  <  h  <  1  ,  and  {|f!uh|||*}0<h<1 
ii)  A:  1/  ->  V  of  (7.6)  is  l/-pseudomonotone 


is  bounded 


(10.5) 

(10.6) 


Proof.  We  follow  the  usual  pseudomonotone  method: 

(1)  a  priori  bounds,  (2)  passage  to  the  limit  and  (3)  the  pseudomono tonicity 


argument. 


1)  From  the  proof  of  the  coercivity  property  of  A  ,  (7.11),  it  is 


apparent  that  A  is  also  coercive  from  V  V'  : 

<A(v),v>_>  aJVv|2  +  Y1 1|  v[|  P  -  Y2  »  v  v  £  v 


(10.7) 


Hence,  by  integrating  equation  (10.1)  with  respect  to  time  from  0  to  T  €  [0,T] 
and  using  formula  (2.9)  and  (10.7),  we  obtain 


\  |Uh(T)|2  -  \  |ujM2  +  Y x  j  ||  Uh(t)  II  P  dt  -  Y2  T 


<  <f(t),Uh(t)>dt 
•'0 


;P'bp’  L1 


f(t)||  P'  dt  +  ^  f  II  Uh(t)  II  P  dt 

v  Jo 


Then,  by  choosing  b  >  0  such  that  Y^  -  bP/p  >  0  ,  it  follows  that  the  sequence 

l  n 

(U  ^Q<h<l  is  bounded  ln  ^  and  in  L°°(0,T;L  (ft) )  . 

2)  With  the  previous  result  and  the  boundedness  of  A  from  1/  -*■  l/'  given 


10.4 


by  (7.8),  the  validity  of  (10.3)  follows  via  weak  compactness  arguments  and, 

then,  upon  the  passage  to  the  limit  in  equation  (10.1),  (10.4)  is  easily  concluded. 

All  of  this  proceeds  as  in  the  case  of  monotone  parabolic  problems;  cf.  [11,  Chap.  2], 

3)  It  remains  to  be  proved  that,  if  either  (10.5)  or  (10.6)  holds,  then 

[X,v]  =  [A(u),v]  ,  V  v  6U  (10.8) 

From  (10.1),  (10.3)  and  (10.4),  we  see  that 

lim  { [Uh,Uh]  +  [A(Uh),Uh]}  =  lim  [f,Uh]  =  [f,u] 
hlO  h+0 

=  [u,u]  +  [X,u] 

=  lim  {[u,Uh]  +  [A(Uh),u]} 
h+0 


Therefore, 

lim  [A(Uh),Uh-u]  =  -  lim  [i^-u.U11]  =  -  j  lim  |  Uh  (T)-u  (T)  | 2 
h!0  h!0  ^  h+0 


5  0  (10.9) 

Now,  by  identical  arguments  to  those  given  in  the  proof  of  (5.5),  (10.8)  follows 
from  (10.9)  and  (10.3)^  when  assuming  either  (10.5),  and  using  the  W-pseudomonotonicity 
property  of  A:  V  ->  I/'  (cf.  Theorem  7.3),  or  (10.6).  This  completes  the  proof 
of  the  theorem. 


We  next  show  that  condition  (10.5)  is  also  sufficient  for  the  strong  convergence 
of  the  approximation  process. 

Theorem  10.3.  Suppose  the  condition  (10.5)  holds  with  bound  y'  >  0  .  Then 
the  subsequence  ^u^o<h<l  Faedo-Galerkin  approximate  solutions  converging 
weakly  to  a  solution  u£ll)  of  problem  (7.5),  in  the  sense  of  Theorem  10.2,  is 
such  that,  as  h  -*■  0+ 


Uh  •+  u  strongly  in  L°°(0,T;L2(f2) ) 
U  -»  u  strongly  in  V 


(10.10) 


In  fact,  the  following  approximation  error  estimates  hold  V  Z  €  LP(0,T;V  ): 

h 


10.5 


|u(T)-Uh(T)|  <  |uo-u"|  +  K1(T,p)|j  u-U 


+  K2(T»Vi»U*  )  I))  u-Z  | 


1/2 


u-U 


hn  p:/2 

LP(Q) 

v  T  e  [0,T] 

hu  l/(p-l) 


(10.11) 


-K3lVUol2/P  +  k4(t,u)||u-u  „ 

J  °  °  *  LP(Q) 


+  K5(T,y,M') HI u-Z | 


|1/P 


(10.12) 


where  p  >  0  is  a  bound  for  u  and  {u  in  1/ 

U<n<l 


Proof .  By  using  formual  (2.9)  and  the  Garding-type  inequality  (7.16)  in 
LP(0,t;W^,P(Q))  ,  t  €.  [ 0 , T]  ,  it  follows  that 

<u(t)-Uh(t)  +  A(u(t) )  -  A(Uh(t)),u(t)-Uh(t)>dt 


>  \  |u(T)-Uh(T)|2  -  \  I uq-Uq 1 2  +  «1 


u(t)-Uh (t)  II  dt 


-  a  00 


u(t)-Uh(t)  ||  P  dt 


LP(fi) 


0 

p'/p 


(10.13) 


and,  from  equations  (7.5)  and  (10.1),  the  following  orthogonality  condition  holds: 


r  t 


<u(t)-U  (t)  +  A(u(t) )  -  A(Uh(t),Z(t)>dt  =  0 


V  z6l/(,T;Vh)  (10.14) 

Hence,  introducing  (10.14)  into  (10.13)  and  using  the  local  Lipschitz  continuity 
property  (7.7),  we  obtain 


\  |u(x)-Uh(T)|2  +  a  [  ||  u(t)-Uh(t)ji  Pdt 

J0 

-  2  iUo~Uoi2  +  a9^')U  “'^ll  PD 


Lf(Q) 


+  <c(p)  1 1! vi— u 


1-^111 


u-Z 


v  x  e  [0.T]  ,  V  Z  e  U(0,T;Vh)  (10.15) 

Therefore,  the  approximation  error  estimates  (10.11)  and  (10.12)  are  implied  by 
(10.15).  Note  that  the  strong  convergence  of  U*1  ■+  u  in  LP(Q)  is  a  consequence 
of  (10. 3) n  assumption  (10.5)  and  the  compact  embedding  of  W  into  LP(Q) 


1 


10.6 


(cf.  Theorem  6.1).  D 

The  Potential  Case.  We  conclude  this  section  by  showing  that  if  the 

O 

bounded,  coercive,  locally  Lipschitz  continuous,  Garding-type  operator  A  of 
(7.6),  is  potential  in  the  sense 

CIII.  A  is  the  gradient  of  some  Gateaux  differentiable  functional  J:  V  R  , 
for  which  there  is  a  constant  y  >  0  such  that 

J(v)  _>  Y ||  v ||  P  ,  V  v  £  V  (10.16) 

then,  for  data 


(f,uo) e  L"(Q)  x  V 


Uh  ->  u 


strongly  in  V 


(10.17) 


(10.18) 


o  o 

the  Faedo-Galerkin  sequence  of  approximations  defined  uniquely  by  (10.1)  is 
such  that 


{Uh}  is  bounded  in  L  (0,T;V) 

U<h<l 

is  bounded  in  L2(Q) 


(10.19) 


Since  L  (Q)C+  V1  ,  property  (10.19),,  is  stronger  than  (10.5)  and,  consequently, 
the  results  of  Theorems  10.2  and  10.3  are  true  in  this  potential  case. 

We  now  prove  this  result  and  establish  the  corresponding  regularity  of 
limit  functions. 

Theorem  10.4.  Let  the  operator  A  of  (7.6)  satisfy  condition  (CIII)  and 
consider  problems  (7.5)  and  (10.1)  with  data  (10.17),  (10.18).  Then  the  Faedo- 
Galerkin  sequence  of  approximate  solutions  ^uh^0<h<l  *S  bounded  in  the  sense 
of  (10.19).  Furthermore,  there  is  a  subsequence  of  approximations,  also  denoted 
{Uh}g<h<1  ,  converging  strongly  to  a  solution  u  of  problem  (7.5)  in  the 

sense  of  (10.10),  such  that,  as  h  ■+  0+  , 


U“ 

& 


weakly*  in  L  (0,T;V) 
weakly  in  L2(Q) 


(10.20) 


10.7 


Proof.  Let  {uh}o<h<i  be  the  Faedo-Gaierkin  sequence  defined  uniquely  by 
(10.1),  (10.17),  (10.18),  which  approximates  problem  (7.5)  with  data  (10.17),  and 

•  Jj 

suppose  that  condition  (CII1)  holds.  Then,  by  replacing  i|;  by  U  in  equation 
(10.1),  integrating  with  respect  to  time  from  0  to  t  €  [0,T]  and,  then, 
observing  that  dJ(Ull(t) )  /dt  =  <£(Uh(t))  ,Uh(t)£>  and  (f(t),U^(t))  £  1/2  ( f  ( t)  J  ^ 

+  1/2 1 Uh (t)  |2  for  a.e.  t  6(0,T),  we  obtain 

if  |i)h(t))2dt  +  y||  uh(T)li  p  <  J(t£)  +  |  !  f  (t)  1 2dt 

z  Jo  Jo 

Vie  [0,T]  (10.21) 


But,  from  the  boundedness  of  A  as  a  map  from  V  V'  (cf.  Remark  7.1), 


><l£)  =  Jo  +  l0<^A(8U^),1^>d8  <  Jo  +  1  11  A(sUo)ll*dSll  ^11  1 


const. 


Therefore,  (10.19)  is  true. 


Next  observe  that  from  Theorem  10.3,  there  is  a  subsequence  of  approximations 
{U^Jq  that  converges  strongly  to  a  solution  u  of  problem  (7.5)  in 

Uf)L  (0,T;L2(£2))  .  Hence,  Uh — *  u  weakly  in  (/  (C^.L1(0,T;V)  densely)  and 
this  together  with  (10.19^  is  equivalent  to  (10.20)^.  Also  (Uh}0<h<1  is 
bounded  in  U  (C*  l /  densely)  and  this  with  U*1 — ^u  weakly  in  V  is  necessary 
and  sufficient  for  Uh — ku  weakly  in  U  (cf.  [17,  Sec.  V.l]).  Then  (10.20)2 
necessarily  holds  and  this  completes  the  proof  of  the  theorem.  D 
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CONCLUSIONS 

For  the  nonlinear  evolution  problems  considered  here,  we  have  shown  that 
coerciveness  and  W-pseudomonotonicity  of  A:  \J  -*■  (/'  guarantees  the  existence 
of  solutions  of  (3.1)  whereas  condition  (8.1)  implies  uniqueness  of  solutions. 

The  elliptic  regularization  ideas  discussed  in  Section. 4  provide  a  general  frame¬ 
work  for  Galerkin  approximations  of  W-pseudomonotone  problems.  We  have  established 
criteria  for  the  existence  and  weak  convergence  of  such  approximations  in  Sections 

O 

4,  5  and  9  and  strong  convergence  whenever  a  Garding-type  inequality  of  the  type 
in  (7.16)  holds.  More  generally,  our  approximation  results  in  Section  9  apply 

O 

to  operators  satisfying  inequalities  of  Garding-type.  In  particular,  the 
existence  and  weak  convergence  of  such  approximations  were  proved  in  Section  4, 

5  ad  9  and,  from  the  analysis  of  Section  9,  it  follows  that  their  convergence 

O 

is  strong  provided  A  satisfies  a  nonlinear  Garding-type  inequality  of  the  form 

[A(v)-A(w)  ,v-w]  _>  a  ||  v-w||  y  -  H(y,  ||  v-w||  ) 

LP(0,T;X) 

V  v ,w  e  (0)  C  V 

where  a  j  >  0  ,  and  X  is  a  Banach  space  continuously  embedded  in  H  and  in 
which  V  is  compactly  embedded.  Also,  if  in  addition.  A:  V  •*  V  is  locally 
Lipschitz  continuous,  we  have  shown  that  error  estimates  for  Galerkin  approximations 
can  be  derived. 

The  Faedo-Galerkin  method  was  considered  as  an  alternative  method  for 
constructing  approximate  solutions.  In  these  cases,  coercivity,  boundedness  and 
demicontinuity  of  A  from  V  -*■  V’  are  sufficient  conditions  for  existence,  and 
local  Lipschitz  continuity  from  V  ■>  V*  is  a  sufficient  condition  for  uniqueness. 

As  we  have  seen,  the  convergence  of  this  method  is  a  conditional  property  in  the 
case  that  A  be  non-monotone.  In  general,  we  may  conclude  the  following  convergence 
results:  Let  A  satisfy  the  existence  conditions  for  the  abstract  problem  and 
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its  Faedo-Galerkin  approximation:  the  conditions  discussed  in  Sections  3  and  10.  Then 
the  Faedo-Ga 1 erk i n  method  is  weakly  convergent  if  (i)  the  sequence  of  time  derivatives 
of  the  approximate  solutions  is  bounded  in  I/'  ,  or  if  (ii)  A:  (/  -*■  V'  is 

pseudomonotone.  The  convergence  of  the  method  is  strong  if  (iii)  condition 

O 

(i)  holds  and  A  satisfies  a  nonlinear  Garding  inequality  of  the  type  given  above. 
Furthermore,  in  the  case  in  which  condition  (iii)  is  satisfied,  error  estimates 
are  derivable  which  are  compatible  with  the  interpolation  theory  of  finite-elements 
in  Sobolev  spaces  [12],  [3]. 

A  fundamental  convergence  condition  for  the  Faedo-Galerkin  method  when 
applied  to  U/-pseudomonotone  parabolic  problems  is  that  the  sequence  time  derivatives 
of  the  approximate  solutions  be  bounded  in  V’  .  We  have  shown  that  this  condition 
is  satisfied  whenever  A  is,  in  addition,  continuous  and  potential  from  V  -+•  V’  , 
its  potential  is  coercive,  and  the  data  (f ,uq) £  H  x  V  .  In  this  potential  case, 

the  convergence  condition  holds  in  H  V  .  Furthermore,  the  approximate  solutions 

00  - 

form  a  sequence  bounded  in  L  (0,T;V)C*V  .  Then  the  regularity  in  time  result 
"(u,9u/3t)  €  L  (0,T;V)  x  H"  holds  for  the  exact  solutions  of  the  problem. 
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APPENDIX  D 


A  Pseudo-Parabolic  Variational  Inequality 


and  Stefan  Problem 


1. 


1 .  Introduction. 

Let  A  and  6  be  maximal  monotone  operators  and  let  C  be  a  non-empty 
closed  convex  set  in  the  real  Hilbert  space  H.  We  shall  give  existence  and 
uniqueness  results  for  evolution  inequalities  (formally)  of  the  form 

(1.1. a)  u(t)  c  C:  (^~(Au(t) )  +  Bu(t)-f(t),  v-u(t))R  >  0  ,  v  e  C  , 

0  <  t  <  T  , 

(1.1. b)  (Au (0) -Vq  ,  v-u (0) ) H  >  0  ,  v  e  C  , 

2 

when  f  e  L  (0,T;H)  and  v  c  Afa^)  are  given.  In  Section  2  we  introduce  a 
new  notion  of  weak  solution  of  (1.1)  and  verify  uniqueness  when  A  is  linear 
self-adjoint  and  B  is  str ictly-monotone .  Existence  of  a  weak  solution  is 
proved  when  A  is  sLrongly-monotone,  B  is  a  subgradient,  and  both  operators 
are  locally  bounded. 

Variational  inequalities  of  die  form  (l . 1 )  are  of  interest  on  their  own  as 
extensions  of  corresponding  evolution  equations  of  Sobolev  type  (where  C =  H) . 

Early  work  on  such  inequalities  is  described  in  [2];  we  mention  [6]  specifically 
as  a  source  of  examples  of  initial-boundary  value  problems  for  the  pseudo-para¬ 
bolic  partial  differential  equation 

(1.2)  — (u-aAu)  =  kAu 

ut 

with  a  >  0,  k  >  0.  Such  equations  arise  as  models  for  diffusion,  and  they  pro¬ 
vide  an  interesting  alternative  to  the  classical  diffusion  equation  wherein  a=0. 
In  Section  3  we  recall  the  two-temperature  heat  conduction  model  from  [3]  and  de¬ 
velop  a  corresponding  one-phase  Stefan  problem  for  (1.2).  Then  we  show  that  such 
a  problem  leads  to  the  variational  inequality  (1.1).  This  development  is  parallel 
to  that  of  the  classical  case  a  =  0  which  is  described,  e.g.,  in  [7].  Existence 
of  a  classical  solution  of  a  Stef. in  problem  for  (1.3)  in  one  dimension  was  given 
in  | 9]  by  entirely  different  methods. 


2. 


2 .  The  Variational  Inequality. 

2 

We  denote  by  L  (0,T;H)  the  Hilbert  space  of  (Bochner)  square-integrable 
functions  on  the  interval  (0,T)  with  values  in  the  Hilbert  space  H.  Let 
1I^(0,T;H)  denote  the  absolutely  continuous  11-valued  functions  v  whose  deriv¬ 
atives  —  belong  to  L2(0,T;H).  Denote  the  dual  of  H  by  H*  and  recall  the 

2  2  * 

natural  identification  L  (0,T;H)  =  L  (0,T;H)  ;  thus  we  obtain  the  (dual)  iden¬ 
tification  L^(0,T;H)c->-  H^(0,T;H)  by  restriction.  The  derivative 
d  1  2 

— -  :  H  (0,T;H)  ->  L  (0,T;H)  is  a  bounded  linear  operator  which  determines  the  dual 
dt 

operator  Ls-(^-)  :  L^(0, T;H)  >  H^(0,T;H)  by  the  formula 


(Lf,v)  = 2 


f  e  L2(0,T;H)  ,  v  c  H1 (0, T ;H) 


The  restriction  of  Lf  to  H-valued  test  functions  is  the  (distribution)  derivative 
~  .  Moreover,  for  f  e  H1(0,T;H)  we  have 


<Lf,v)  =  (: 


djf 

dt 


v)  9 


(f(0),v(0))u-  (f(T),v(T))H 


v  e  Hl(0,T;H) 


Thus,  we  can  regard  "Lf +f(T)"  as  formally  equivalent  to  the  Cauchy  operator 
"  7T  +  f(0)  "  . 


We  shall  use  basic  material  on  maximal  monotone  operators  (1].  Specifically, 

recall  A  c  HxH  is  monotone  if  (x^  ,  y^]  e  A  for  j  =  l  and  2  imply 

(Xi-X2  ,  ^1^2^  —  strictly  monotone  if  in  addition  equality  holds  only  if 

x.  = x„  ,  and  strongly  monotone  if  there  is  a  c  >  0  for  which  (x.-x  ,  y,-y,)„  > 
2 

c||x1-x2H^  for  all  such  pairs  [x^  ,  y  ].  If  ip:  H  *  R  U  {  +<*>)  is  proper,  convex 
and  lower  semicontinuous,  its  sub grad  lent  defined  by 


()<p(x)  ?  (u  c  H:  (u, y-x)  <  tp(y)-ip(x)  for  all  y  e  H) 

H 


3. 


for  x  £  H  is  maximal  monotone.  More  specifically,  il  C  is  a  non-empLy,  convex 
and  closed  set  in  H,  its  indicator  i unc t ion 

{0  ,  x  e  C 

+  ,  x  i  C 

is  proper,  convex  and  lower  semicontinuous,  and  we  have  u  e  dl  (x)  if  and  only 

v/ 

if 

x  £  C:  (u,y-x)  <  0  for  all  y  £  C  . 

We  shall  be  concerned  with  maximal  monotone  operators  A  with  domain  D(A)  *  H. 

That  is,  A(x)  is  non-empty  for  every  x  e  H.  This  is  known  to  be  equivalent 

to  A  being  locally  bounded :  for  each  x  e  H  there  is  a  neighborhood  of  x  on 

which  A  is  bounded.  This  does  not  imply  A  is  bounded  in  general  unless  H 

has  finite  dimension  or  A  is  linear. 

Suppose  we  are  given  a  pair  A,B  of  maximal  monotone  operators  on  the 

2 

Hilbert  space  H,  a  closed  convex  subset  C  of  H,  f  c  L  (0,T;H)  and 
[Uq  ,  Vg]  e  A.  The  triple  (u,v,w)  is  a  strong  solution  of  (1.1)  if 

u, v  e  H^O.T-.H);  w  e  L2(0,T;H)  , 

(2.1. a)  u(t)  e  C:  (~^*  +  w  (t) -f  (t ) ,  x-u(t)>H  >  0  ,  x  e  C  , 

v(t)  e  A(u(t))  and  w(t)  c  B(u(t))  for  a.e.  t  e  [0,T]  ,  and 

(2 . 1  .b)  (v(0)-vQ  ,  x-u(0))H  >  0  ,  x  e  C  . 

Mote  that  since  u  and  v  are  continuous,  C  is  closed  in  H  and  A  is 
closed  in  HxH,  it  follows  that  the  inclusions  u(t)  e  C  and  v(t)  e  A(u(t)) 
hold  for  all  t  e  [0,T] .  Also,  (2.1)  can  be  restated  as 


4. 


(2. 2. a)  +  w(t)  +  OIc(u(t))  3  f(t) 

(2 .2 .b)  v(0)  +  <.Ic(u(0))  3  vQ 

in  terms  of  the  indicator  function. 

We  shall  relax  the  requirement  that  v  e  H^(0, T;H)  as  follows.  Set 
K  £  {u  e  H* (0,T;H) :  u(t)  e  C,  0  <  t  <  T) .  Define  a  weak  solution  of  (1 .1)  to  be 
a  triple  {u,v,w}  satisfying 

u  e  K  ;  v,w  e  L2(0,T;H)  , 

v (t)  e  A(u(t) )  ,  w(t)  e  B  (u (t) )  ,  a.e.  t  e  [0,T]  , 
and  for  some  ?  e  A(u(T))  we  have 

(2.3)  (Lv+w-f.q-u)  +  (?,,’)  (T)-u(T))h  >  (vQ  ,  q(0)-u(0))H  ,  rj  e  K  . 

Note  that  if  {u,v,w]  is  a  strong  solution  then  it  is  a  weak  solution  with 
£=v(T).  Moreover  we  have  the  following  elementary  result. 

Theorem  1.  Let  A  be  continuous,  linear,  self-adjoint  and  monotone;  let  B  be 
strictly  monotone.  Then  the  first  two  components  of  a  weak  solution  are  uniquely 
determined. 


Proof :  Let  (u^  ,  ,  v^)  be  weak  solutions  for  j  =  1,2.  By  our  assumptions  on 

A  we  may  assume  (after  modification  on  a  null  set)  v  ^  =  AUj  o  H^(0,T;H)  and 
^=A(u^(T)).  Thus  we  have 

<LAut  +w1-f,u2-u1>  +  (Au1(T),u2(T)-ul(T))H  >  (vQ  ,  u2 (0)-^ (0)>H 
(LAu2+w2-f,ul-u2>  +  (Au2(T),u1(T)-u2(T))H  >  (vQ  ,  ^  (0)  (0)  )R  . 


<LAu,u>  =  1/2((Au(0),u(0))h-(Au(T),u(T))h)  , 

so  adding  the  two  inequalities  and  applying  this  identity  with  u*u^-u2  gives 

(w.-w,  ,  u.-u,)  „  +  1/2 (Au (T),u(T))  +  1/2  (Au(0)  ,u(0))„  <  0  . 

1212  L2(0,T;H)  H  H 

Strict  monotoneity  of  B  shows  u^  =  u2  • 

Remark •  Without  additional  assumptions  on  the  set  C  we  should  not  expect  any 

uniqueness  of  the  third  component,  w.  For  example,  in  the  extreme  case  C  =  {0), 

2 

(2.3)  is  vacuous  and  we  need  only  choose  v,w  e  L  (0,T;H)  with  v(t)  e  A(0)  and 

w(t)  e  B  (0)  to  obtain  a  weak  solution.  On  the  other  hand,  if  C  =  H  then  any 

dv  2 

weak  solution  is  a  strong  solution  of  the  equation  — +w=f  In  ^  (0,T;H)  with 
initial  condition  v(0)=Vq  .  See  [5]  for  such  Cauchy  problems. 

The  primary  objective  here  is  the  following  existence  result. 

Theorem  2.  Let  C  be  a  non-empty,  closed  convex  set  in  Hilbert  space  H.  Let 

A  and  B  be  maximal  monotone  operators  on  H  such  that  A  is  strongly  monotone 

with  domain  D(A)=H,  B  is  a  subgradient,  B  =  '<p,  with  D(B)  * H  and  <p(x)  ^0 

2 

for  all  x  e  H.  For  u^  e  C,  v^  e  A(Uq),  and  f  c  L  (0, T;H)  given,  there  is  at 
least  one  weak  solution  (u,v,w)  of  (1.1). 

Remarks .  Since  A  is  strongly  monotone  we  may  assume  it  is  of  the  form  A  +  I. 

Thus  we  wish  to  replace  (2.3)  by 

(2.4)  <L(u+v)4w-f,T1-u)  +  (u(T)+?,'](T)-u(T))H  >  (u0  +  vQ  ,  i(0)-u(0))H  ,  n  e  K  . 

Proof :  For  each  c  0  let  be  the  Yoshida  approximation  of  the  indicator 

function  I  .  The  subdifferential  i'I*  is  a  maximal  monotone  Lipschitz  continuous 

w  w 


6. 


function  on  H  and  we  have  n(tj>  +  l£)  =o<p+ol£  .  From  Theorem  1  of  [5]  we  obtain 
a  strong  solution  of  the  approximating  Cauchy  problem 

(2.5)  ^r(ue(t)  +vc(t))  +  VE(t)  +  01p(uc(t))  =  f(t>  , 

(t)  e  A(uE(t))  ,  wc<t)  e  B(u£(t))  ,  a.e.  t  e  [0,T]  , 

ut(0)  =  uo  *  V0)  "  v0  e  A<uo)  * 

1  2 

with  u  ,  v  e  H  (0,T;H)  and  w  e  L  (0,T;H).  Taking  the  inner  product  with 

C  C  c 

du 

*T—  in  (2.5)  and  using  the  chain  rule  | 1,  Lemma  3-3]  give 


I 


t  du 

I  — 

.  1  dt 


dr  + 


I 


(We  dropped  the  non-negative  term 


t  dv  du 

I(df  •ir)H " 


by  monotoneity  of  A.) 
du  « 


Therefore 


du 


—If  +  <p(u  (t))  +  lUu  (t))  <  [|f||  ll-rfll  , 

L  (0, T ;H)  C  C  C  h  (0,T;H)  dt  LZ(0,T;H) 


+  «p(uQ)  +I^(u0) 


Recall  that 


(2.6)  ijj(x)  =  ^-|x-Projc(x)  |Z  ,  x  e  H  ; 

thus  I^(ut(t))  >  0  and  I^(uQ)-0  since  uQ  e  C.  From  <j>(x)  >  0,  x  e  H,  we 
obtain 


sup  |u  (t)  |  + 

0<t<T 


(0,  T ;H) 


(2.7) 


<  M 


(2.8) 


sup  I(u(t))<M 
0<t<T  L  E 


where  M  depends  on  and  f  but  not  e.  By  the  Ascoli-Arzela  Theorem  we  may 

pass  to  a  suitable  subnet  (indexed  again  by  e)  and  obtain 


(2-9) 


(2.10) 


u^(t)  *  u(t)  ,  strongly  in  H  ,  uniformly  in  t  , 

du  d  2 

•JjT  +  57  .  weakly  in  L  (0,T;H)  . 


The  limit  u  e  H  (0,T;H)  is  continuous  so  its  range  is  a  compact  path  in  H. 
Since  A  and  B  are  locally  bounded  there  is  an  e ^  such  that  for  0  <  c  < 
and  t  e  [0,TJ 

Ive  (t)  lH  +  fwt(t>lH  • 

Thus  we  may  use  the  maximality  of  A  and  B  to  pass  to  a  subnet  (again  indexed 
by  c)  for  which 


(2.11) 


V£(T)  S  weakly  in  H  ,  and 


(2.12)  ■>  v  ,  wt  *  w  weakly  in  LT(0,T;H)  , 


where  5  e  A(u(T))  and  v(t)  c  A(u(t)),  w(t)  <  B(u(t))  for  almost  every  t  e  [0,TJ 
We  shall  show  the  triple  {u,v,w]  is  a  weak  solution.  From  (2.6),  (2.8) 


follows 


|uc(t)-Projc(uE(t))  |H  <  2eM  ,  t  e  [0,T)  , 
so  (2.9),  (2.10)  show  u  e  K.  For  any  rj  e  K  we  obtain  from  (2.5) 

<l(u£  +vc)  +we-f,  +  (Uf;  (T)  +  Ve  (T),  tj  (T)-ue  (T))r 

=  (-dIE(u  M-u)  +  (u  (0) +v  (0),rj(0)-u  (0))„  . 

^  e  C  l/(0,T;H)  t  c  c  H 
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From  the  definition  of  subgradient  and  the  inclusion  t)  e  K  we  obtain 


(-01! (u  ),.i-u ) ,  >  r  dj(u  (t))-i^(T,(t))dt  >  0 , 

C  C  e  L  (0, T;H)  J0  C  E  C 


so  there  follows 

<L(u£ +Vc) +we-f,r1-uc)  +  (ue(T) +V£(T),t1(T)-ue(T))h  >  (uQ  +v0  ,  ")  (°)  ~ue  (°)  )H  • 

Using  (2. 9) -(2. 12)  and  the  weak  continuity  of  L,  we  take  the  limit  as  e  •>  0 
in  the  preceding  inequality  and  obtain  (2.4). 


Remark.  The  approximation  of  (2.4)  by  (2.5)  is  an  abstract  penalty  method. 


3.  A  Stefan  Problem. 

We  consider  a  problem  of  heat  diffusion  involving  a  solid-liquid  phase 
change  at  a  prescribed  temperature.  One  application  we  have  in  mind  is  the 
melting  of  ice  (initially  at  temperature  zero)  suspended  in  a  reservoir  or  porous 
medium.  The  novelty  in  this  treatment  is  that  we  assume  the  heat  diffusion  is 
governed  by  the  pair  of  equations 

~  =  kAip  ,  9  =  9-  aAip  . 

Chen  and  Gurtin  [3]  introduced  such  a  model  for  heat  conduction  in  non-simple 
materials  where  the  energy,  entropy,  heat  flux  and  thermodynamic  temperature 
8(x, t)  depend  on  the  conductive  temperature  <p(x, t)  and  its  first  two  spatial 
gradients.  Here  the  heat  flux  is  determined  by  the  conductive  temperature  and 
the  phase  is  determined  by  the  thermodynamic  temperature.  Thus  9  >  0  in  the 
region  occupied  by  water  and  9=0  corresponds  to  the  frozen  region. 


9. 


We  describe  the  geometry  of  the  problem.  Let  the  bounded  domain  G  in  Rn 
be  the  medium  in  which  the  ice/water  is  suspended  and  let  its  boundary  dG  consist 
of  two  disjoint  pieces,  Tq  and  r  .  Set  ll«Gx  (0,T),  where  T  >  0,  and  note 

that  its  lateral  boundary  is  Bq  u  B^  ,  where  B^  =  P^  x  (0,  T)  for  J*0, 1.  The 
water-region  ft^  =  {(x,t)  e  ft:  9(x,t)  >  0}  is  separated  from  the  Ice-region 

{  (x,  t)  e  ft:  9(x,t)=0)  by  an  interface  S  which  is  the  phase  boundary.  The 
unit  outward  normal  on  is  denoted  by  N=  (N^  ,  N^),  e  R°.  If  V(t)  is 

the  velocity  in  Rn  of  the  interface  at  time  t,  then  it  follows  by  the  chain 

rule  that  V(t)-Nx+Nt  =  0  on  S.  Set  n  =  / IlN^ II >  the  unit  outward  normal  in 

n  -*» 

R  of  the  lateral  boundary  of  .  Of  course  n=Nx  on  »  and  =*  0  where 

t  =  0  or  t  =  T. 

The  problem  is  formulated  as  follows.  We  are  given  the  conductivity  k  >  0, 
temperature  discrepancy  a  >  0,  and  latent  heat  b  >  0,  of  the  material  and  a 
constant  h  >  0  representing  conductivity  across  the  lateral  boundary  .  The 
initial  thermodynamic  temperature  9Q(x),  x  e  G,  and  applied  conductive  temper¬ 
ature  g(x,t),  (x,t)  e  ,  are  given  with  9Q=0  on  VQ  ,  9Q  >  0  on  r  ,  and 

g  >  0.  The  local  form  of  the  problem  is  to  find  a  pair  of  non-negative  functions 

9,«p  on  i.  for  which  we  have 


(3-D 

^  =  kA*p  ,  9  =  4>  -  aAip 

in 

U1 

(3-2) 

k  +  bV(t)-n  =  0 

on 

on 

S 

(3.3) 

k  In  +  =  0 

on 

ri 

(3.4) 

ip  =  0 

on 

9( ' , 0)  =  9q 


(3.5) 


on  G  . 


•  *  ' 


r* 
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Note  that  if  0,  ip  is  a  solution  of  (3.1)-(3.5)  and  0^  >  0,  then 

ciQ 

(3.6)  (a/k)  ^  +  0  =  v  in  w 

so  it  follows  that  <p=0  on  ;.q  u  S.  Since  g  >  0,  the  maximum  principle  for 
the  elliptic  equation  in  (3.1)  on  the  region  G(t)  s{x  e  G:  (x,  t)  e  shows 

that  (p  >  0  in  and  <  0  on  S.  Thus  Nt  <  0  on  S  and  G(t)  is  in¬ 

creasing  with  t. 

We  shall  show  that  the  problem  (3.1)-(3.5)  leads  to  a  variational  inequality 

of  the  form  (1.1).  Letting  H^(G)  denote  the  Sobolev  space  of  functions  v  in 

L^(G)  for  which  all  derivatives  ,  1  <  j  <  n,  belong  to  L^(G),  we  define 

i  j 

V  =  (v  c  H  (G) :  v|  =0}.  Here  v|  is  the  trace  on  the  boundary  of  G;  see 

ro  0 

[8,10)  for  details.  Regarding  regularity  of  a  solution,  we  assume  0^  e  V, 

0:  [0, T]  *  V  is  absolutely  continuous,  ip  e  L*(0,T;V),  and  (c.f.  (3.6)) 

(3*7)  £•  +  8(t)  -  <p(t)  ,  a . e.  t  e  [0,T]  . 

★ 

Define  the  continuous  linear  B:  V  *  V  by 


Bu 


(v)  =  k(Vu-  v"v)dx  +  J*  h(uv)ds  ,  u,v  e  V 


For  a  test  function  v  e  Cq((0,T),V)  we  obtain 


T 

J*  Btp(t)  (v(t)  )dt=J  k  V«p-  W  dxdt  +  J  h  <p  v  dsdt 


J*  (-kAip)v  dxdt  +  J*  k  Vtp-N^vdsdt  +  J*  htpvdsdt 

"J  It V  +  ?  h8v  +  [  bNtv 

dt  JB,  JS  f 


from  (3.1)- (3.4).  Furthermore  we  have 


J.vl 


oH(e) 

nt 


(V) 


in  the  sense  of  V  -valued  distributions,  where  H(s)  *  1  for  s  >  0  and 
H(s)  =  0  for  s  <  0  is  the  Heaviside  function.  We  can  summarize  the  above 
calculations  as 


(3.8)  -^-(8+bH(e))  +  Bo  *  (hg),,  in  Ll(0,T;V*) 

at  i  j 

where  we  define 


v  e  V  , 


t  e  [  0,  T]  . 


Combining  (3.7)  and  (3.8)  we  find  Chat  the  absolutely  continuous  function 
9:  [ 0, T]  *  V  satisfies 

(3. 9. a)  ~(0  +  (a/k)B (0)  +bH(8) )  +  B(9)  =  (gh)p  in  L^O.TjV*)  , 
dt  I\ 

(3. 9-b)  0(0)  =  90  ,  and 

(3.9.c)  9(x,t)  >  0  ,  a.e.  x  e.  C  ,  t  e  [0,T]  . 


If  we  integrate  (3. 9. a)  and  set 


u(t)  f  0(s)ds 
J0 


f(t)  (I  +  (a/k)B  +  bll)9  -  b  +  f  (hg)p  (s)ds  , 

0  J0  ‘l 


there  follows 
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(3.10)  “7(1 +  (a/k)B)u  +  Bu  -  f (t)  =  b(l-H(®))  . 

at 

Finally  we  note  that  H(u)  =  H(®)  since  G(t)  is  increasing  in  t,  hence, 
u(l-H(®) )  =  0  in  w. 

The  preceding  computations  show  that  u  e  H^(0,T;V)  and  it  satisfies 
u(0)  =  0, 

fu(t)  >0  in  V  , 

^■(1+  (a/k)B)u(t)  +Bu(t)  >  f(t)  in  V*.  and 
(~(I  +  (a/k)B)u(t)  +Bu(t)-f (t))(u(t))  =  0,  0<t<T  . 

Setting  C  =  (v  e  V:  v  >  0  a.e.  in  G}  we  see  that  u  is  a  strong  solution  of 

(1.1)  with  A=I  +  (a/k)B  and  Uq=v0=0;  c.f.  (2.1).  Note  that  we  can  trivially 

★ 

rephrase  the  material  of  Section  2  in  the  H  -  H  duality  instead  of  the  H  -  H 
pairing  through  the  scalar  product. 

Theorem  1  asserts  uniqueness  of  a  solution  of  (3.1)- (3. 5)  under  conditions 
considerably  weaker  than  those  leading  to  (3.11).  Theorem  2  establishes  existence 
of  a  weak  solution  to  (3.11)  which  possesses  certain  additional  regularity  prop¬ 
erties.  These  topics  will  be  developed  elsewhere  by  other  methods  [4,11]. 


13. 
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